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Introduction 
Topics Covered in this Chapter 


e Prerequisites and text books 

e Scalar, vector and tensor fields 

e Curves, surfaces, and volumes 

e Coordinate systems 

e Units 

e Continuum approximation 

e Densities, potential gradients, and fluxes 
e Velocity: a measure of flux by convection 
e Density 

e Species concentration 

e Energy (heat) 

e Porous media 

e Momentum 

e Electricity and Magnetism 


Introduction 


This course is designed as a first level graduate course in transport 
phenomena. Undergraduate courses generally start with simple example 
problems and lead to more complex problems. With this approach, the student 
must learn the fundamental principles by induction. The approach used here 
is to teach the fundamental principles and then deduce the analysis for 
example problems. The example problems are classical problems that should 
be familiar to all Ph.D. Chemical Engineering graduates. These problems will 
be presented not only as an exercise with analytical or numerical solutions but 
also as simulated experiments which are to be interpreted and graphically 
displayed for presentation. 


Prerequisites and text books 


Students in this class are expected to have a background corresponding to a 
BS degree in Chemical Engineering. This includes a course in multivariable 
calculus, which covers the algebra and calculus of vectors fields on volumes, 


surfaces, and curves of 3-D space and time. Courses in ordinary and partial 
differential equations are a prerequisite. Some elementary understanding of 
fluid mechanics is expected from a course in transport phenomena, fluid 
mechanics, or physics. It is assumed that not all students have the prerequisite 
background. Thus, material such as vector algebra and calculus will be briefly 
reviewed and exercise problems assigned that will require more reading from 
the student if they are not already familiar with the material. 


The two required textbooks for this course are R. Aris, Vectors, Tensors, 
andthe Basic Equations of Fluid Mechanics and Bird, Stewart, and Lightfoot, 
Transport Phenomena. Several of the classical problems are from S. W. 
Churchill, Viscous Flows, The Practical Use of Theory. The classical 
textbook, Feyman, Leighton, and Sands, The Feyman Lectures on Physics, 
Volume II is highly recommended for its clarity of presentation of vector 
fields and physical phenomena. The students are expected to be competent in 
MATLAB, FORTRAN, and EXCEL and have access to Numerical Recipes in 
FORTRAN. 


The following table is a suggested book list for independent studies in 
transport phenomena. 


Author Title Publisher Year 

L.D. Landau , ia 

andE. M. Fluid Mechanics, 2 Butterworth 1987 
peste Ed. 

Lifshitz 

V. G. Levich Physicochemical Prentice-Hall 1962 

Hydrodynamics 
S Hydrodynamics and 
Chaativacdenge Hydromagnetic Dover 1961 


Stability 


Author 


H. Schlichting 


H. Lamb 


S. Goldstein 


W. E. Langlois 


J. Happel, H. 
Brenner 


G. K. 
Batchelor 


S.-I. Pai 


M. Van Dyke 


S. W. 
Churchill 


S. W. 
Churchill 


R. F. Probstein 


Title 


Boundary Layer 
Theory 


Hydrodynamics 
Modern 
Developments in 
Fluid Dynamics 
Slow Viscous Flow 
Low Reynolds 
Number 


Hydrodynamics 


An Introduction to 
Fluid Mechanics 


Viscous Flow 
Theory I Laminar 
Flow 
Perturbation 


Methods in Fluid 
Mechanics 


Inertial Flows 


Viscous Flows 


Physicochemical 
Hydrodynamics 


Publisher 


McGraw-Hill 


Dover 


Dover 


Macmillan 


Kluwer 


Cambridge 


Van Nostrand 


Academic 
Press 


Etaner 


Butterworths 


Butterworth- 
Heinemann 


Year 


1960 


1932 


1965 


1964 


1973 


1967 


1956 


1964 


1980 


1988 


1989 


Author 


S. Middleman 


S. Middleman 


E. L. 
Koschmeider 


W.-J. Yang 
W.-J. Yang 


A. J. Chorin 


Le. 
Wrobel,C. A. 
Brebbia 


M. J. 
Baines,K. W. 
Morton 


Title 


An Introduction to 
Fluid Dynamics 


An Introduction to 
Mass and Heat 
Transfer 


Benard Cells and 
Taylor Vortices 


Handbook of Flow 
Visualization 


Computer-Assisted 
Flow Visualization 


Computational Fluid 
Mechanics 


A Mathematical 
Introduction to Fluid 
Mechanics 


Computational 
Modeling of Free 
and Moving 
Boundary Problems, 
Vol. 1 Fluid Flow 


Numerical Methods 
for Fluid Dynamics 


Publisher 


John Wiley 


John Wiley 


Cambridge 


Taylor & 
Francis 


CRC Press 


Academic 
Press 


Springer- 
Verlag 


Computational 
Mechanics 
Publications 


Oxford 


Year 


1998 


1998 


1993 


1989 


1994 


1989 


1993 


1991 


1993 


Author 


W. E. 
Schiesser,C. 
A. Silebi 


N. Ida, J. P. A. 


Bastos 


L. G. Leal 


W. M. Deen 


C.S. Jog 


C.S. Jog 


T.J. Chung 


R. J. Kee, 


M.E. Coltrin, 


P. Glarborg 


Title 


Computational 
Transport 
Phenomena 


Electro-Magnetics 
and Calculation of 
Fields 


Laminar Flow and 
Convective 
Transport Processes 


Analysis of 
Transport 
Phenomena 


Foundations and 
Applications of 
Mechanics, Vol. I, 
Continuum 
Mechanics 


Foundations and 
Applications of 
Mechanics, Vol. II, 
Fluid Mechanics 


Computational Fluid 
Dynamics 


Chemically Reacting 
Flow 


Publisher 


Cambridge 


Springer 


Butterworth 


Oxford 


CRC Press 


CRC Press 


Cambridge 


Wiley - 
Interscience 


Year 


1997 


1997 


1992 


1998 


2002 


2002 


2002 


2003 


Author Title Publisher Year 


Fluid Dynamics; 


7.U.A. Warsi Theoretical and Taylor & 2006 
Computational Francis 
Approaches 
Y. A. Cengel Fluid Mechanics; 
and J. M. Fundamentals and McGraw Hill 2006 
Cimbala Applications 


Transport phenomena book list 


Scalar, vector and tensor fields 


Scalars, vectors, and matrices are concepts that may have been introduced to 
the student in a course in linear algebra. Here, scalar, vector, and tensor fields 
are entities that are defined over some region of 3-D space and time. It is 
implicit that they are a function of the spatial coordinates and time, i.e., 

y = »(2, y, z,t) = y(a, t). The spatial coordinates are expressed as 
Cartesian coordinates in this class. However, vectors and tensors are physical 
entities that are independent of the choice of spatial coordinates even though 
their components depend on the choice of coordinates. 


Scalar fields have a single number, a scalar, at each point in space. An 
example is the temperature of a body. The temperature field is usually 
expressed visually by a contour map showing curves of constant temperature 
or isotherms. An alternative visual display of a scalar field is a color map 
with the value of the scalar scaled to a gray scale, hue, or saturation. The 
values of the scalar field are continuous with the exception of definable 
surfaces of discontinuity. An example is the density of two fluids separated 
by an interface. Media that are chaotic and discontinuous on a microscopic 
scale may be described by an average value in a representative elementary 
volume that is large compared to the microscopic heterogeneity but small 
compared to macroscopic variations. An example is the porosity of a porous 
solid. 


Vector fields have a magnitude and direction associated with each point in 
space. An example is the velocity field of a fluid in motion. Vector fields in 
two dimensions can be visually expressed as field lines that are everywhere 
tangent to the vector field and whose separation quantifies the magnitude of 
the field. Streamlines are the field lines of the velocity field. Alternatively, a 
vector field in two dimensions can be visually expressed by arrows whose 
directions are parallel to the vector and having a width and/or length that 
scales to the magnitude of the vector. These graphical representations of 
vector fields are not useful in three dimensions. In general, a vector field in 3- 
D can be expressed in terms of its components projected on to the axis of a 
coordinate system. Thus, a vector field may have different components when 
projected on to different coordinate systems. Since a vector is a physical 
entity, the components in different frames of reference transform by 
prescribed rules. The position of a point in space relative to an origin is a 
vector defined by the distance and direction. Special vectors having a 
magnitude of unity are called unit vectors and are used to define a direction 
such as coordinate directions or the normal direction to a surface. We will 
denote vectors with bold face letters, e.g., v, x, or n. 


Tensors are physical entities associated with two directions. For example, the 
stress tensor represents the force per unit area, each of which are directional 
quantities. The velocity gradient is a tensor. Transport coefficients, such as 
the thermal conductivity, are tensors, which transform a potential gradient to 
a flux, each of which are vectors. The components of a tensor in a particular 
coordinate system are represented by a 3x3 matrix. Since the tensor is a 
physical entity that is independent of the coordinate system, the components 
must satisfy certain transformation rules between coordinate systems. In 
particular, a set of three directions called the principal directions can be found 
to transform the components of the tensor to a diagonal matrix. This 
corresponds to finding the eigenvectors of a matrix and the components 
correspond to the eigenvalues. Bold face letters will also denote tensors. The 
stress tensor will be denoted by T or t, depending on whether discussing Aris 
or BSL, respectively. 


Curves, surfaces, and volumes 


We will be dealing with regions of space, V, having volume that may be 
bound by surfaces, S, having area. Regions of the surface may be bound by a 
closed curve, C, having length. 


Surfaces are defined by one relationship between the spatial coordinates. 
Equation: 


g? = fiean")s GPE (are a") = 0, or F(x) =0 


Alternatively, a pair of surface coordinates, u!, uw? can define a surface. 
Equation: 


zg’ = 2° (u', uw’), 1=1,2.3.07° x =x (a aie) 


Each point on the surface that has continuous first derivatives has associated 
with it the normal vector, n, a unit vector that is perpendicular or normal to 
the surface and is outwardly directed if it is a closed surface. Fluid-fluid 
interfaces need to also be characterized by the mean curvature, H, at each 
point on the surface to describe the normal component of the momentum 
balance across the interface. The flux of a vector, f, across a differential 
element of the surface is denoted as follows, i.e. the normal component of the 
flux vector multiplied by the differential area. 

Equation: 


f-nda 


Curves are defined by two relationships between the spatial coordinates or by 
the intersection of two surfaces. 
Equation: 


fi (aa as) = Oand f, (eae) = 0, or f, (x) = O0and f, (x) = 0 
Alternatively, a curve in space can be parameterized by a single parameter, 


such as the distance along the curve, s or time, t. 
Equation: 


x= x(5) 


The tangent vector is a unit vector that is tangent to each point on the curve. 
Equation: 


t = dx(s)/ds 


The component of a vector, f, tangent to a differential element of a curve is 
denoted as follows. 
Equation: 


f-tds 


If the parameter along the curve is time, the differential of position with 
respect to time is the velocity vector and the differential of velocity is 
acceleration. 

Equation: 


Coordinate systems 


Scalars, vectors, and tensors are physical entities that are independent of the 
choice of coordinate systems. However, the components of vectors and 
tensors depend on the choice of coordinate systems. The algebra and calculus 
of vectors and tensors will be illustrated here with Cartesian coordinate 
systems but these operations are valid with any coordinate system. The 
student is suggested to read Aris to learn about curvilinear coordinate 
systems. Bird, Stewart, and Lightfoot express the components of the relevant 
vector and tensor equations in Cartesian, cylindrical polar, and spherical polar 
coordinate systems. 


Cartesian coordinates have coordinate axes that have the same direction in the 
entire space and the coordinate values have the units of length. Curvilinear 
coordinates, in general, may have coordinate axis that are in different 
directions at different locations in space and have coordinate values that may 
not have the units of length, e.g., 6 in the cylindrical polar system. If (y!, y?, 
y’) are Cartesian coordinates and (z', x’, x?) are curvilinear coordinates, a 
differential length is related to the differential of the coordinates by the 
following relations. 


3 
ds? = > dy" dy* 
k=1 


dy* = ov dx' = ov da! + oy da? ae OF da’ 
3 k k 
Oy ; Oy 
ds? = -dx' ——dzx/ 
va On! OxI 


where gj are components of the metric tensor which transforms differential of 
the coordinates to differential of length. Summation is understood for 
repeated indices. Calculus in a curvilinear coordinate system will require the 
metric tensor. 


A differential element of volume in curvilinear coordinate system is related to 
differentials of the coordinates by the square root of the determinate of the 
metric tensor or the Jacobian, J. 
Equation: 
dV = dy' dy’ dy? 
= OY. 204) “09%, of 198 
oe eS dx! dx* dx3 
= g/? dx! dx? dx? 
= Jdzx' dx? dx* 


Henceforth, Cartesian coordinates with subscript notation will be used. 


Units 


Dimensional quantities will be used in equations without explicit 
specification of units because it is understood that they will have the SI 
system of units. The SI units and mks units are similar with some exceptions 
as in electricity and magnetism. The following table lists the SI units of the 
quantities used in this course and the conversion factor needed to convert the 
quantity from some customary units to SI units. Multiply the quantity in 
customary units by the conversion factor to obtain the quantity in SI units. 
The following is taken from, The SI Metric System of Units and SPE METRIC 
STANDARD), Society of Petroleum Engineers. 


Customary Conversion 
Quantity SI unit unit factor 
Length m ft 3.048 E-01 
Mass kg Ibm 4.535 924 E- 

01 
Time S S 1.0 
Temperature °K °R 5/9 
; 6.894 757 

Pressure, stress Pa psi E+03 


Density kg/m? g/cm? 1.0 E+03 


Customary Conversion 

Quantity SI unit unit factor 
4.448 222 

Force N Ibf E+00 
Biguuete may US. . 6.309 020 E- 

gal/min 05 
Diffusivity m?/s cm?/s 1.0 E-04 
Thermal Btu/(hr-ft?- 1.730 735 
conductivity WH) | cE E+00 
Heat transfer 2, Btu/(hr-ft- 5.678 263 E- 
coefficient Wert) °F) 06 
Permeability m? darcy aa a 
Surface tension N/m dyne/cm 1.0 E-03 
Viscosity . : 
(dynamic) Paes cp 1.0 E-03 


SI units and conversion factors 


Continuum approximation 


The calculus of scalar, vector, and tensor fields require that these quantities be 
piecewise continuous down to infinitesimal dimensions. However, quantities 
such as density, pressure, and velocity become ambiguous or stochastic at the 
scale of molecular dimensions. Thus the fields discussed here are the average 
value of the quantity over a representative elementary volume, REV, of space 
that is large compared to molecular dimensions but small compared to the 
macroscopic variation of the quantities. The size of the REV depends on the 


scale that a problem is being investigated. For example, suppose one is 
investigating a fixed bed catalyst reactor. As a first order approximation for 
design purposes, the reactor may be modeled as a one-dimensional system 
with the cross-section of the reactor approximated as the REV. However, is 
one is investigating instabilities and channeling, the bed may be modeled in 
2-D with the REV being a volume that is small compared to the macroscopic 
dimensions of the reactor but large compared to the size of the catalyst 
particles. If one is optimizing the transport-limited kinetics of the reactor, 
then the REV may be small compared to the size of the catalyst particle. If 
one is optimizing the balance between transport-limited and surface reaction 
rate limited kinetics, the REV may be small enough to describe the surface 
morphology of the catalyst particle. However, the molecular dynamics of the 
surface reaction is beyond the realm of transport phenomena. 


Densities, potential gradients, and fluxes 


Velocity: and flux by convection. Transport or flux of the various quantities 
discussed in this course will be due to convection (or advection) or due to the 
gradient of a potential. Common to all of these transport process is the 
convective transport resulting from the net or average motion of the 
molecules or the velocity field, v. The convective flux of a quantity is equal 
to the product of the density of that quantity and the velocity. In this sense, 
the velocity vector can be interpreted as a “volumetric flux” as it has the units 
of the flow of volume across a unit area of surface per unit of time. Because 
the flux by convection is common to all forms of transport, the integral and 
differential calculus that follow the convective motion of the fluid will be 
defined. These will be known as the Reynolds’ transport theorem and the 
convective or material derivative. 


Mass density and mass flux. If p is the mass density, the mass flux is pv. 


Species concentration. Suppose the concentration of species A in a mixture is 
denoted by Cy. The convective flux of species A is Cy v. Fick’s law of 
diffusion gives the diffusive flux of A. 

Equation: 


J,=-D,eVC, 


The diffusivity, Da, is in general a tensor but in an isotropic medium, it is 
usually expressed as a scalar. 


Internal energy (heat). The density of internal energy is the product of density 
and specific internal energy, pE. The convective flux is pE v. For an 
incompressible fluid, the convective flux becomes p Cp (T-To) v. The 
conductive heat flux, q, is given by Fourier’s law for conduction of heat, 
Equation: 


q= —-ke VT 


where k is the thermal conductivity tensor (note: same symbol as for 
permeability). 


Porous media. The density of a single fluid phase per unit bulk volume of 
porous media is @p, where @ is the porosity. Darcy’s law gives the volumetric 
flux, superficial velocity, or Darcy’s velocity as a function of a potential 
gradient. 

Equation: 


u = —7 ¢(Vp— pg) 
= ov 


where k is the permeability tensor and v is the interstitial velocity or the 
average velocity of the fluid in the pore space. Darcy’s law is the momentum 
balance for a fluid in porous media at low Reynolds number. 


Momentum balance. Newton’s law of motion for an element of fluid is 
described by Cauchy’s equation of motion. 
Equation: 


pa=pa 
=pf{+VeT 


where f is the sum of body forces and T is the stress tensor. The stress tensor 
can be interpreted as the flux of force acting on the bounding surface of an 
element of fluid. 

Equation: 


| 
<J 
e 
4 
Q 
SN 


I (oa — pf)dv 
V(t) V(t) 


| 
= 

| 

5 

Q 
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The stress tensor for a Newtonian fluid is as follows. 
Equation: 


where p is the thermodynamic pressure, © is the divergence of velocity, pL is 
the coefficient of shear viscosity, (A+2/3y) coefficient of bulk viscosity, and e 
is the rate of deformation tensor. Thus the anisotropic part (not identical in all 
directions) of the stress tensor is proportional to the symmetric part of the 
velocity gradient tensor and the constant of proportionality is the shear 
viscosity. 


Electricity and Magnetism. We will not be solving problems in electricity and 
magnetism but the fundamental equations are presented here to illustrate the 
similarity between the field theory of transport phenomena and the classical 
field theory of electricity and magnetism. The Maxwell’s equations and the 
constitutive equations are as follows. 

Equation: 


VeD=po 
Constitutive equations: 
B= uH 

D = cE 

J =oE 


where 


¢ E electric field intensity 

D electric flux density or electric induction 

e H magnetic field intensity 

¢ B magnetic flux density or magnetic induction 
e J electric current density 

e p charge density 

¢ 1 magnetic permeability (tensor if anisotropic) 
e eelectric permittivity (tensor if anisotropic) 

e oelectric conductivity (tensor if anisotropic) 


When the fields are quasi-static, the coupling between the electric and 
magnetic fields simplify and the fields can be represented by potentials. 
Equation: 


EK = —-VV 
B=VxA 


where V is the electric potential and A is the vector potential. The electric 
potential is analogous to the flow potential for invicid, irrotational flow and 
the vector potential is analogous to the stream function in two-dimensional, 
incompressible flow. 


Reading assignment 


Note:Read Chapter 1 and Appendix A and B of Aris. 


Cartesian Vectors and Tensors: Their Algebra 
Topics Covered in this Chapter 


e Definition of a vector 

e Examples of vectors 

e Scalar multiplication 

e Addition of vectors — coplanar vectors 
e Unit vectors 

e A basis of non-coplanar vectors 

e Scalar product — orthogonality 

e Directional cosines for coordinate transformation 
e Vector product 

e Velocity due to rigid body rotations 

e Triple scalar product 

e Triple vector product 

e Second order tensors 

e Examples of second order tensors 

e Scalar multiplication and addition 

e Contraction and multiplication 

e The vector of an antisymmetric tensor 
e Canonical form of a symmetric tensor 


Note:Reading Assignment: Chapter 2 of Aris, Appendix A of BSL 


The algebra of vectors and tensors will be described here with Cartesian 
coordinates so the student can see the operations in terms of its components 
without the complexity of curvilinear coordinate systems. 


Definition of a vector 


Suppose 2; , i.e., (%1, £2, £3), are the Cartesian coordinates of a point P in 
a frame of reference, 0123. Let 0123 be another Cartesian frame of 
reference with the same origin but defined by a rigid rotation. The 


coordinates of the point P in the new frame of reference is x; where the 
coordinates are related to those in the old frame as follows. 
Equation: 


= 45 4 11524 fe lo; a 13523 


ej, = 1; 2; = 121 + let, + ligx3 


where 1; are the cosine of the angle between the old and new coordinate 
systems. Summation over repeated indices is understood when a term or a 
product appears with a common index. 


Cartesian Vector 
A Cartesian vector, a, in three dimensions is a quantity with three 
components @1, @2, a3 in the frame of reference 0123, which, under 
rotation of the coordinate frame to 0123 , become components 
Q1,@2,a3 , where 


Equation: 


aj= 150; 


Examples of vectors 


In Cartesian coordinates, the length of the position vector of a point from 
the origin is equal to the square root of the sum of the square of the 
coordinates. The magnitude of a vector, a, is defined as follows. 
Equation: 


lal = (aja;)”? 


A vector with a magnitude of unity is called a unit vector. The vector, a/|al, 
is a unit vector with the direction of a. Its components are equal to the 
cosine of the angle between a and the coordinate axis. Some special unit 


vectors are the unit vectors in the direction of the coordinate axis and the 
normal vector of a surface. 


Scalar multiplication 


If a is a scalar and a is a vector, the product aa is a vector with components, 
Qj, Magnitude aja|, and the same direction as a. 


Addition of vectors — Coplanar vectors 


If a and b are vectors with components a; and b;, then the sum of a and b is 
a vector with components, a;+b;. 


The order and association of the addition of vectors are immaterial. 
Equation: 


at+b=b+a 
(a+b)+c=a+(b+c) 


The subtraction of one vector from another is the same as multiplying one 
by the scalar (-1) and adding the resulting vectors. 


If a and b are two vectors from the same origin, they are colinear or parallel 
if one is a linear combination of the other, i.e., they both have the same 
direction. If a and b are two vectors from the same origin, then all linear 
combination of a and b are in the same plane as a and D, i.,e., they are 
coplanar. We will prove this statement when we get to the triple scalar 
product. 


Unit vectors 


The unit vectors in the direction of a set of mutually orthogonal coordinate 
axis are defined as follows. 


Equation: 


The suffixes to e are enclosed in parentheses to show that they do not 
denote components. A vector, a, can be expressed in terms of its 
components, (a,, d>,a3) and the unit vectors. 

Equation: 


a = a1€(1) + A2€(2) + A3€(3) 
This equation can be multiplied and divided by the magnitude of a to 


express the vector in terms of its magnitude and direction. 
Equation: 


a = |al (Ze + Te) + fee) 
= |al (Area) + Aze(a) + Aseg)) 


where A; are the directional cosines of a. 


A special unit vector we will use often is the normal vector to a surface, n. 
The components of the normal vector are the directional cosines of the 
normal direction to the surface. 


Scalar product — Orthogonality 


The scalar product (or dot product) of two vectors, a and b is defined as 
Equation: 


ae b = |a||b| cosé 


where 8 is the angle between the two vectors. If the two vectors are 
perpendicular to each other, i.e., they are orthogonal, then the scalar 
product is zero. The unit vectors along the Cartesian coordinate axis are 
orthogonal and their scalar product is equal to the Kronecker delta. 
Equation: 


eeey = di 
flings 
— ((0,t 4G 


The scalar product is commutative and distributive. The cosine of the angle 
measured from a to b is the same as measured from b to a. Thus the scalar 
product can be expressed in terms of the components of the vectors. 
Equation: 


aeb = (a1e(1) ae a2&(2) + a3€(3)) e (bye(1) a b2€(2) = b3e(3)) 


ab ji; 


| 


a,b; 


| 


The scalar product of a vector with itself is the square of the magnitude of 
the vector. 


Equation: 
aea = |a||alcos0 
= |al’ 
aea = a4;a; 
= |al’ 


The most common application of the scalar product is the projection or 
component of a vector in the direction of another vector. For example, 
suppose n is a unit vector (e.g., the normal to a surface) the component of a 
in the direction of n is as follows. 

Equation: 


aen = |a|cos6 


Directional Cosines for Coordinate Transformation 


The properties of the directional cosines for the rotation of the Cartesian 
coordinate reference frame can now be easily illustrated. Suppose the unit 
vectors in the original system is e,;) and in the rotated system is €,,) . The 


components of the unit vector, €(;) , in the original reference frame is ;;. 


This can be expressed as the scalar product. 
Equation: 


eg) = Lijec + laje(ay + I3je(s),5 = 1,2,3 
E(i) © &Gj) = ligt _ ie 2 3 


Since €(,) is a unit vector, it has a magnitude of unity. 
Equation: 


egy eg) = 1 =liyligy = Ly + layla) + layla I = 1, 2,3 


Also, the axis of a Cartesian system are orthorgonal. 


Equation: 
0,ift #7 
€(i) ° &G) = we =j 
thus 
ei) © eG) = dij 
Equation: 


e() © ej) = leiley = Lislij + leila; + lsilsj, 7,7 = 1, 2,3 
=e 


Vector Product 


The vector product (or cross product) of two vectors, a and b, denoted as 
axb, is a vector that is perpendicular to the plane of a and b such that a, b, 
and axb form a right-handed system. (i.e., a, b, and axb have the 
orientation of the thumb, first finger, and third finger of the right hand.) It 
has the following magnitude where @ is the angle between a and b. 
Equation: 


la x b| = |al |b] sind 


The magnitude of the vector product is equal to the area of a parallelogram 
two of whose sides are the vectors a and b. 


Since the vector product forms a right handed system, the product bxa has 
the same magnitude but opposite direction as axb, i.e., the vector product is 
not commutative, 

Equation: 


bxa=-axb 


The vector product of a vector with itself or with a parallel vector is zero or 
the null vector, i.e., axa=0. A quantity that is the negative of itself is zero. 
Also, the angle between parallel vectors is zero and thus the sine is zero. 


Consider the vector product of the unit vectors. They are all of unit length 
and mutually orthogonal so their vector products will be unit vectors. 
Remembering the right-handed rule, we therefore have 

Equation: 


The components of the vector product can be expressed in terms of the 
components of a and b and applying the above relations between the unit 
vectors. 

Equation: 


as Db <= (a1e(1) + a2€ (2) + a3€(3)) x (b1e(1) + be (2) + b3e(3)) 
= (azb3 _ a3b2)e(1) ++ (a3b1 — a1b3)e,2) se (a1b2 _ a2b)e,3) 


The permutations of the indices and signs in the expression for the vector 
product may be difficult to remember. Notice that the expression is the 
same as that for the expansion of a determinate of the matrix, 

Equation: 


ay “S@) 2G) 
at ag a3 
b, by 63 


Expansion of determinants are aided by the permutation symbol, €;;x. 
Equation: 


0, if any two of i,7,k are the same 
Eijk = +1, if 17k isan even permutation of 1,2,3 
—1, if 77k is an odd permutation of 1,2,3 


The expression for the vector product is now as follows. 
Equation: 


axhb=— EijkAiD;€(K) 


Velocity due to rigid body rotations 


We will show that the velocity field of a rigid body can be described by two 
vectors, a translation velocity, vj, and an angular velocity, w. A rigid body 
has the constraint that the distance between two points in the body does not 
change with time. The translation velocity is the velocity of a fixed point, 
O, in the body, e.g., the center of mass. Now consider a new reference 
frame (coordinate system) with the origin at point O that is translating with 
respect to the original reference frame with the velocity vq. The rotation of 
the body about O is defined by the angular velocity, o, i.e., with a 
magnitude w and a direction of the axis of rotation, n, such that the positive 
direction is the direction that a right handed screw advances when subject to 
the rotation, .o= n. Consider a point P not on the axis of rotation, having 
coordinates x in the new reference frame. The velocity of P in the new 
reference frame has a magnitude equal to the product of @ and the radius of 
the point P from the axis of rotation. This radius is equal to the magnitude 
of x and the sine of the angle between x and a, i.e., |x| sin®. The velocity of 
point P in the new reference frame can be expressed as 

Equation: 


vVv=WXxxX 


|v| = w |x| sin 8 


The velocity field of any point of the rigid body in the original reference 
frame is now 
Equation: 


V=Vi) tWxx 


Since this equation is valid for any pair of points in the rigid body, the 
relative velocity Av between a pair of points separated by Ax can be 
expressed as follows. 

Equation: 


Av = wx Ax 


Conversely, if the relative velocity between any pair of points is described 
by the above equation with the same value of angular velocity, then the 
motion is due to a rigid body rotation. 


Triple scalar product 


The triple scalar product is the scalar product of the first vector with the 
vector product of the other two vectors. It is denoted as (abc) or [abc]. 
Equation: 


(abc) = ae(b xc) 


Recall that bxc has a magnitude equal to the area of a parallelogram with 
sides b and c and a direction normal to the plane of b and c. The scalar 
product of this normal vector and the vector a is equal to the altitude of the 
parallelepiped with a common origin and sides a, b, and c. The triple scalar 
product has a magnitude equal to the volume of a parallelepiped with a 
common origin and sides a, b, and c. The sign of the triple scalar product 
can be either positive or negative. If a, b, and c are coplanar, then the 
altitude of the parallelepiped is zero and thus the triple scalar product is 
zero. 


The triple scalar product can be expressed in terms of the components by 
using the earlier definitions of the vector product and scalar product. 
Equation: 


oe EijkdiCj€(K) 
a = Am€(m) 
ae (b x c) = EijkAmDiCj€(m) e E(k) = EijkA@mdiCjOmk = EjjkAKDiC; 
= €;jhAibjCK 


From the definition of the permutation symbol, the triple scalar product is 
unchanged by even permutations of a, b, and c but have the opposite 
algebraic sign for odd permutations. Also, if any two of a, b, and c are 
identical, then permutation of the two identical vectors results in a triple 


scalar products that are identical and also opposite in sign. This implies that 
the triple scalar product is zero if two of the vectors are identical. 


Triple vector product 


The triple vector product of vectors a, b, and c results from the repeated 
application of the vector product, i.e., ax(bxc). Since bxc is normal to the 
plane of b and c and ax(bXxc) is normal to bxc, ax(bxc) must be in the 
plane of b and c. It is left as an exercise to show that 

Equation: 


a x (bx c) = (aec)b— (aeb)c 


Second order tensors 


A second order tensor can be written as a 3X3 matrix. 
Equation: 


Ai, Ai Aig 
A= Axi Azo Ags 
Az, Az. A33 


A tensor is a physical entity that is the same quantity in different coordinate 
systems. Thus a second order tensor is defined as an entity whose 
components transform on rotation of the Cartesian frame of reference as 
follows. 

Equation: 


A pq = lipl jqAij 


If Aij=Ajithe tensor is said to be symmetric and a symmetric tensor has only 
six distinct components. If Aij=-Aijthe tensor is said to be antisymmetric 
and such a tensor is characterized by only three nonzero components for the 
diagonal terms, Ai;, are zero. The tensor whose ij element is Aji is called 


the transpose A’ of A. The determinant of a tensor is the determinant of the 
matrix of its components. 
Equation: 


det A = €4;~A1;A2;A3x 


Examples of second order tensors 


A second order tensor we have already encountered is the Kronecker delta 
6;;. Of its nine components, the six off-diagonal components vanish and the 
three diagonal components are equal to unity. It transforms as a tensor upon 
transforming its components to a rotated frame of reference. 

Equation: 


je) 


pq — lip! jq9ij 


Ll; 


aq 


ip 


) 


Pq 


because of the orthogonality relation between the directional cosines J;. In 
fact, the components of 6;in all coordinates remain the same. 6; is called 
the isotropic tensor for that reason. The transport coefficients (e.g., thermal 
conductivity) of an isotropic medium can be expressed as a scalar quantity 
multiplying 6j. 


If a and b are two vectors, the set of nine products, a ; b ; =A jj, is a second 
order tensor, for 


Equation: 


ip 


Ang = apd, = Lipail jqbj = lipl jq (aib;) 
Lint jgAij 


tp" q 


An important example of this is the momentum flux tensor. If p is the 
density and v is the velocity, p v; is the i‘* component in the direction Oi. 
The rate at which this momentum crosses a unit area normal to Oj is the 
tensor, P Vj; Vj. 


Scalar multiplication and addition 


If a is a scalar and A a second order tensor, the scalar product of a and A is 
a tensor aA each of whose components is a times the corresponding 
component of A. 


The sum of two second order tensors is a second order tensor each of whose 
components is the sum of the corresponding components of the two tensors. 
Thus the ij" component of A+B is Ajj + Bj. Notice that the tensors must be 
of the same order to be added; a vector can not be added to a second order 
tensor. A linear combination of tensors results from using both scalar 
multiplication and addition. aA + &B is the tensor whose ij" component is 
aA;; + &B;;. Subtraction may therefore be defined by putting a = 1, 6 = —1 


Any second order tensor can be represented as the sum of a symmetric part 
and an antisymmetric part. For 
Equation: 


1 1 
Aig = 3 (Aig + Agi) + 5 (Ais — Asi) 


and changing i and j in the first factor leaves it unchanged but changes the 
sign of the second. Thus, 
Equation: 


1 1 
A= Z(AtAN + S(A- A) 


represents A as the sum of a symmetric tensor and antisymmetric tensor. 


Contraction and multiplication 


As in vector operations, summation over repeated indices is understood 
with tensor operations. The operation of identifying two indices of a tensor 
and so summing on them is known as contraction. A ;; is the only 
contraction of A 
Equation: 


ip 
Aj, = Ai + Age + Ass 


and this is no longer a tensor of the second order but a scalar, or a tensor of 
order zero. The scalar A;; is known as the trace of the second order tensor A. 
The notation tr A is sometimes used. The contraction operation in 
computing the trace of a tensor A is analogous to the operation in the 
calculation the magnitude of a vector a, i.e., 

jal? =a-a—aja;,+ a2Qa9 + Aazaz 


If A and B are two second order tensors, we can form 81 numbers from the 
products of the 9 components of each. The full set of these products is a 
fourth order tensor. Contracted products result in second order or zero order 
tensors. We will not have an occasion to use products of tensors in our 
course. 


The product Aij aj of a tensor A and a vector a is a vector whose i" 


component is Aj; a;. Another possible product of these two is Aj a; . These 
may be written A:a and a‘A, respectively. For example, the diffusive flux 
of a quantity is computed as the contracted product of the transport 
coefficient tensor and the potential gradient vector, e.g.,.q = —k- VT 


The vector of an antisymmetric tensor 


We showed earlier that a second order tensor can be represented as the sum 
of a symmetric part and an antisymmetric part. Also, an antisymmetric 
tensor is characterized by three numbers. We will later show that the 
antisymmetric part of the velocity gradient tensor represents the local 
rotation of the fluid or body. Here, we will develop the relation between the 


angular velocity vector, @, introduced earlier and the corresponding 
antisymmetric tensor. 


Recall that the relative velocity between a pair of points in a rigid body was 
described as follows. 
Equation: 


Av =wx Ax 
We wish to define a tensor @ that also can determine the relative velocity. 
Equation: 


Av = wx Ax 
= Axel) 


The following relation between the components satisfies this relation. 
Equation: 


25; = C4jkWk 


1 
Wh = 7 Eijk llij 


Written in matrix notation these are as follows. 


Equation: 
Wy 0 W3 —W9 
W= Wo ,= —-w3 ODO Wy 
W3 W9 Wy 0 


The notation vec @ is sometimes used for @. In summary, an antisymmetric 
tensor is completely characterized by the vector, vec @. 


Canonical form of a symmetric tensor 


We showed earlier that any second order tensor can be represented as a sum 
of a symmetric part and an antisymmetric part. The symmetric part is 
determined by 6 numbers. We now seek the properties of the symmetric 
part. A theorem in linear algebra states that a symmetric matrix with real 
elements can be transformed by its eigenvectors to a diagonal matrix with 
real elements corresponding the eigenvalues. (see Appendix A of Aris.) If 
the eigenvalues are distinct, then the eigenvector directions are orthogonal. 
The eigenvectors determine a coordinate system such that the contracted 
product of the tensor with unit vectors along the coordinate axis is a parallel 
vector with a magnitude equal to the corresponding eigenvalue. The surface 
described by the contracted product of all unit vectors in this transformed 
coordinate system is an ellipsoid with axes corresponding to the coordinate 
directions. 


The eigenvalues and the scalar invariants of a second order tensor can be 
determined from the characteristic equation. 
Equation: 


det (Aj; — Ad;;) = py — AG+ 70 — 23 

where 

O = Ay + Ag + Az3 = trA 

® = Ao. A33 — A3A32 + A33A11 — A3iAi3 + A11A22 — Az Aa 
w = det A 


Assignment 2.1 


a. Relative velocity of points in a rigid body. If x and y are two points 
inside a rigid body that is translating and rotating, determine the 
relation between the relative velocity of these two points as a function 
of their relative positions. If x and y are points on a line parallel to the 
axis of rotation, what is their relative velocity? If x and y are points on 
opposite sides of the axis of rotation but with equal distance, r, what is 
their relative velocity? Draw diagrams. 

b. Prove that: ae(bxc) = (axb)ec 

c. Show ae(bxc) vanishes identically if two of the three vectors are 
proportional of one another. 


d. Show that if a is coplanar with b and c, then a*(bxc) is zero. 

e. Prove that: a x (b x c) = (aec)b— (aeb)c 

f. Prove that the contracted product of a tensor A and a vector a, A:a, 
transforms under a rotation of coordinates as a vector. 


g. Show that you get the same result for relative velocity whether you use 
@ or Q for the rotation of a rigid body. 


Cartesian Vectors and Tensors: Their Calculus 
Topics Covered in this Chapter 


e Tensor functions of time-like variable 

e Curves in space 

e Line integrals 

e Surface integrals 

e Volume integrals 

e Change of variables with multiple integrals 
e Vector fields 

e The vector operator -gradient of a scalar 

e The divergence of a vector field 

e The curl of a vector field 

e Green's theorem and some of its variants 

e Stokes' theorem 

e The classification and representation of vector fields 
e Irrotational vector fields 

e Solenoidal vector fields 

¢ Helmholtz’ representation 

e Vector and scalar potential 


Reading assignment: Chapter 3 of Aris 


Tensor functions of time-like variable 


In the last chapter, vectors and tensors were defined as quantities with 
components that transform in a certain way with rotation of coordinates. 
Suppose now that these quantities are a function of time. The derivatives of 
these quantities with time will transform in the same way and thus are tensors 
of the same order. The most important derivatives are the velocity and 
acceleration. 

Equation: 


dx; 


f dx, 
a(t) = x(t),a; = 


The differentiation of products of tensors proceeds according to the usual rules 
of differentiation of products. In particular, 
Equation: 


(aeb) = B@eb+ac® 


a0 db 
(axb)= 4 xbt+ax 


t 


Jn Sa 


Q 


t 


Curves in space 


The trajectory of a particle moving is space defines a curve that can be defined 
with time as parameter along the curve. A curve in space is also defined by the 
intersection of two surfaces, but points along the curve are not associated with 
time. We will show that a natural parameter for both curves is the distance 
along the curve. 


The variable position vector x(t) describes the motion of a particle. For a finite 
interval of t, say a < t < b, we can plot the position as a curve in space. If the 
curve does not cross itself (i.e., if x (t) A x(t'),a <t <t' < D) itis called 
simple; if x(a) = x(b) the curve is closed. The variable t is now just a 
parameter along the curve that may be thought of as the time in motion of the 
particle. If t and t’ are the parameters of two points, the cord joining them is 
the vector x (t’) — x (t). As t > t’ this vector approaches (t' — t)x (t) and so 
in the limit is proportional to x (t). However the limit of the cord is the tangent 
so that x (¢) is in the direction of the tangent. If v? = x (t) e x (t) we can 
construct a unit tangent vector T. 

Equation: 


T= xX) /0=Vv/0 
Now we will parameterize a curve with distance along the curve rather than 


time. If x(t) and x(t + dt) are two very close points, 
Equation: 


x (t+ dt) = x(t) + dtx(t) + O (dt’) 


and the distance between them is 
Equation: 


ds’ = {x(t + dt) — x(t)} e {x(t + dt) — x(t)} 
x (t) e x(t) dt? + O (dt?) 


The arc length from any given point t = a is therefore 
Equation: 


s(t) = / [x (t’) x (t')] "ae! 


s is the natural parameter to use on the curve, and we observe that 
Equation: 


A curve for which a length can be so calculated is called rectifiable. From this 
point on we will regard s as the parameter, identifying ¢ with s and letting the 
dot denote differentiation with respect to s. Thus 

Equation: 


t(s) = x(s) 


is the unit tangent vector. Let x(s), x(s + ds), and x(s — ds) be three nearby 
points on the curve. A plane that passes through these three points is defined by 
the linear combinations of the cord vectors joining the points. This plane 
containing the points must also contain the vectors 
Equation: 

Rlstds) xis) = x(s) + O (ds) 

and 


x(o+ds)—2x(o)}xls-d9) _ ¥ (5) + O (ds?) 


Thus, in the limit when the points are coincident, the plane reaches a limiting 
position defined by the first two derivatives x (s) and X (s). This limiting 
plane is called the osculating plane and the curve appears to lie in this plane in 
the intermediate neighborhood of the point. To prove this statement: (1) A 
plane is defined by the two vectors, x (s) and X (s), if they are not co-linear. 
(2) The coordinates of the three points on the curve in the previous two 
equations are a linear combination of x (s), x (s) and X (s), thus they line in 
the plane. 


Now x = Tso X = Tandsince Te T = l, 
Equation: 


d(tet) 


Pe —O0=Tettirtet—2te07t 


Tet—0 


so that the vector T is at right angles to the tangent. Let 1/ denote the 


magnitude of T. 
Equation: 


Then v is a unit normal and defines the direction of the so-called principle 
normal to the curve. 


To interpret p , we observe that the small angle d@ between the tangents at s 
and s + ds is given by 
Equation: 
cos dO = T(s) e 7(s + ds) 
1— 5d0°+...=rTer+retds+ STetds*+... 


=1-jretds*+... 


since 7e¢7 =OandsotTe7+7¢e7 = 0. Thus, 
Equation: 


is the reciprocal of the rate of change of the angle of the tangent with arc 
length, i.e., g is the radius of curvature. Its reciprocal 1/p is the curvature, 
Kk = |d6/ds| = 1/p. 


A second normal to the curve may be taken to form a right-hand system with 7 
and v . This is called the unit binormal, 
Equation: 


p=TXxXV 


Line integrals 


If F(a) is a function of position and C is a curve composed of connected arcs 
of simple curves, x = x(t),a <t < borx = x(s),a < s < b, we can define 
the integral of F along C’ as 

Equation: 


So F(x)dt = i. F|x 


or 
fo F(x)ds = J? Fix(o){& (0) ee ()} "at 
Henceforth, we will assume that the curve has been parameterized with respect 


to distance along the curve, s. 


The integral is from a to b. If the integral is in the opposite direction with 
opposite limits, then the integral will have the same magnitude but opposite 


sign. If x(a) = x(b), the curve C is closed and the integral is sometimes 
written 
Equation: 


f Flx(s)las 


If the integral around any simple closed curve vanishes, then the value of the 
integral from any pair of points a and 6 is independent of path. To see this we 
take any two paths between a and 6, say C; and C3, and denote by C' the 
closed path formed by following C’; from a to 6 and C’2 back from 6 to a. 
Equation: 


fo F ds 


[[fFas|] + [Lh aslo, 
[Je Fas]] , — [seas] | 


= 0 


C2 


If a(x) is any vector function of position, a e 7 is the projection of a tangent to 
the curve. The integral of a e 7 around a simple closed curve C is called the 
circulation of a around C. 

Equation: 


f aetds = f ai [2 (8), 22 (8), 25 (s)] 7% ds 


We will show later that a vector whose circulation around any simple closed 
curve vanishes is the gradient of a scalar. 


Surface integrals 


Many types of surfaces and considered in transport phenomena. Most often the 
surfaces are the boundaries of volumetric region of space where boundary 
conditions are specified. The surfaces could also be internal boundaries where 


the material properties change between two media. Finally the surface itself 
may the subject of interest, e.g. the statics and dynamics of soap films. 


A proper mathematical treatment of surfaces requires some definitions. A 
closed surface is one which lies within a bounded region of space and has an 
inside and outside. If the normal to the surface varies continuously over a part 
of the surface, that part is called smooth. The surface may be made up of a 
number of subregions, which are smooth and are called piece-wise smooth. A 
closed curve on a surface, which can be continuously shrunk to a point, is 
called reducible. If all closed curves on a surface are reducible, the surface is 
called simply connected. The sphere is simply connected but a torus is not. 


If a surface is not closed, it normally has a space curve as its boundary, as for 
example a hemisphere with the equator as boundary. It has two sides if it is 
impossible to go from a point on one side to the other along a continuous curve 
that does not cross the boundary curve. The surface is sometimes called the cap 
of the space curve. 


If S is a piece-wise smooth surface with two sides in three-dimensional space, 
we can divide it up into a large number of small surface regions such that the 
dimensions of the regions go to zero as the number of regions go to infinity. If 
the regions fill the surface and are not overlapping, then sum of the areas of the 
regions is equal to the area of the surface. If the function, F' is defined on the 
surface, it can be evaluated for some point of each subregion of the surface and 
the sum ’F' AS computed. The limit as the number of regions go to infinity 
and the dimensions of the regions go to zero is surface integral of F' over S. 
Equation: 


lim °F AS = [ [Fas 
S 


The traditional symbol of the double integral is retained because if the surface 
is a plane or the surface is projected on to a plane, then Cartesian coordinates 
can be defined such that the surface integral is a double integral of the two 
coordinates in the plane. Also, two surface coordinates can define a surface and 
the double integration is over the surface integrals. 


In transport phenomena the surface integral usually represents the flow or flux 
of a quantity across the surface and the function F’ is the normal component of 
a vector or the contracted product of a tensor with the unit normal vector. Thus 
one needs to know the direction of the normal in addition to the differential 
area to calculate the surface integral. Consider the case of a surface defined as 
a function of two surface coordinates. 

Equation: 


To see how we arrive at this result, recall the partial derivatives of the 
coordinates of a curve with respect to a parameter is a vector that is tangent to 
the curve. The magnitude is 

Equation: 


Of — of e of 7 
Ou; ~~ Ou; Ou; 


The vector product has a magnitude equal to the product of the magnitudes and 
the sine of the angle between the vectors. This gives us the area of a 
parallelogram corresponding to the area of the differential region. 


Equation: 
as = (z=) (z=) sin 0 du, dug 
du, u2 duy ral 


The two tangent vectors in the direction of the surface coordinates lie in the 
tangent plane of the surface. Thus the direction of the vector product is 


perpendicular to the surface. Inward or outward direction for the normal has 
not yet been specified and will be determined by the sign. 


Volume integrals 


The volume integral of a function F’ over a volumetric region of space V is the 
limit of the sum of the products of the volume of small volumetric subregions 
of V and the function F' evaluated somewhere within each subregion. 


Equation: 
[[ [Fe dV =lim °F AV 
V 


Change of variables with multiple integrals 


In Cartesian coordinates the elements of volume dV is simply the volume of a 
rectangular parallelepiped of sides dx, dx2, dx3 and so 
Equation: 


dV = dx, dx2 dx3 


Suppose, however, that it is convenient to describe the position by some other 
coordinates, say £1, €2, £3. We may ask what volume is to be associated with 
the three small changes d£1, dé2, d&3. 


The change of coordinates must be given by specifying the Cartesian point x 
that is to correspond to a given set €1, £9, €3, by 
Equation: 


Li = Lj (1, £9, &3) 


Then by partial differentiation the small differences corresponding to a change 
d€; are 


Equation: 


= Ox; 


dz; = 
0§; 


d§ j 


Let dx") be the vectors with the components (6x; /d€;)d&; for 7 = 1,2, and3. 
Then the volume element is 


Equation: 
dV = dxe (dx) » dx'3)) 
Lj Oz ; x 
Cijk oe BE, 182 ees 
= Jd&; d&2 d&3 

where 

Equation: 

O(x1, £2, £3) Ox, Ox; Ox, 


a O(&1, £9, £3) me 06, by Obs, 


is called the Jacobian of the transformation of variables 


Vector fields 


When the components of a vector or tensor depend on the coordinates we 
speak of a vector or tensor field. The flow of a fluid is a perfect realization of a 
vector field for at each point in the region of flow we have a velocity vector 
v(x). If the flow is unsteady then the velocity depends on the time as well as 
position, v = v(x, t). 


Associated with any vector field a(x) are its trajectories, which is the name 
given to the family of curves everywhere tangent to the local vector a. They 
are solutions of the simultaneous equations 

Equation: 


d dx; 
OF a (x); thatis—* = a; (x1, 2, r3) 
ds ds 


Where s is a parameter along the trajectory. (It will be arc length if a is always 
a unit vector.) Streamlines of a steady flow are a realization of these 
trajectories. For a time dependent vector field the trajectories will also be time 
dependent. If C’ is any closed curve in the vector field and we take the 
trajectories through all points of C,, they describe a surface known as the vector 
tube of the field. For flow fields, it is called a stream tube. 


The vector operator V-gradient of a scalar 


The symbol V (enunciated as "del") is used for the symbolic vector operator 
whose i” component is 5/dz;. Thus if V operates on a scalar function of 
position $(x) it produces a vector Vd with components 6¢/62;. 

Equation: 


Ob Ob On; ad 


Oz ; Ox; On; 7 0 Oa; 


Equation: 


6) ra) Fa) 
poe eee ay = e) ea +e) a 


since a; = 1;;a; so V@¢ is a vector. The suffix notation, 2 for the partial 
derivative with respect to x; is a very convenient one and will be taken over for 
the generalization of this operation that must be made for non-Cartesian frame 
of reference. The notation "grad" for V is often used and referred to as the 
gradient operator. Thus grad ¢ is the vector which is the gradient of the scalar. 
The gradient operator can also operate on higher order tensors and the 
operation raises the order by one. Thus the gradient of a vector a is grad a, Va, 
or in component notation a;,;. V is sometimes written 5/dx and can be 
expanded as the vector operator 


The suffix notation, z for the partial derivative with respect to x; is a very 
convenient one and will be taken over for the generalization of this operation 
that must be made for non-Cartesian frame of reference. The notation "grad" 
for V is often used and referred to as the gradient operator. Thus grad ¢ is the 
vector which is the gradient of the scalar. The gradient operator can also 
operate on higher order tensors and the operation raises the order by one. Thus 
the gradient of a vector a is grad a, Va, or in component notation a;,;. V is 
sometimes written 6 /6x and can be expanded as the vector operator 
Equation: 


V=e. 0 re) te O de 0 
eae |) ae. Ox; bo) ar Ax1 (2) “a Axs (3) 5 Ox 


since a; = 1;;a; so V@ is a vector. 


Equation: 
0 
dé = 32 da; 
a a a 
= on dx, = a ig dx3 


= (ea ge +e +eG $2) © (eq der + e@ dea + e() des) 


The unit vector in the direction of dx is u = dx/ds. The derivative of ¢ in the 
direction of u is 
Equation: 

dp 


— = V¢ eu = |V¢| cos 0 
ds 


If (a) = cis a surface, then V¢ is normal to the surface. To prove this, let dx 
be a differential distance on the surface. The differential of ¢ along dx is zero 
for any dx on the surface. This implies that the scalar product of V@ with any 
vector on the surface is zero or that V@ has zero component or projection on 
the tangent plane and thus V¢ is normal to the surface. Also, since V¢ is 


normal to the surface, the derivative of ¢ is a maximum in the direction normal 
to the surface. 


The divergence of a vector field 


The symbolic scalar or dot product of the operator V and a vector is called the 
divergence of the vector field. Thus for any differentiable vector field a(x) we 
write 
Equation: 

jaa VeeHogs 3 
Or 0x2 0x3 


The divergence is a scalar because it is the scalar product and because it is the 
contraction of the second order tensor a; ;. 


We will now demonstrate why the Ve operation on a vector field is called the 
divergence. Suppose that an elementary parallelepiped is set up with one corner 
P at x1, £2, £3 and the diagonally opposite one Q at x; + dx 1, r2 + dx, 

x3 + dx3 as shown in Fig. 3.7. The outward unit normal to the face through Q 
which is perpendicular to 01 is e(;), whereas the outward normal to the parallel 


face through P is —e(). Suppose a(x) is a differentiable flux vector field. We 


are going to compute the net flux of a across the bounding surfaces of the 
parallelepiped. The value of the normal component of a at some point on the 
two faces perpendicular to the 01 direction are 

Equation: 


al (x14, L2; x3)anda, (x1 + dx1, £25 £3) 
where 
Ko < ©) < 42 + drgandx3 < 73 < 43 + dx3 


Thus if n denotes the outward normal and dS is the area dx2dz3 of these two 
faces, we have a contribution from them to the surface integrak fan dS of 


Ss 
Equation: 
Sf aendS = [ay (x1 + dx}, L225 £3) — ay (x1, L225 r3)| dxo dx3 
Olfaces 


— oa dx, dzydz3 +O (dx* ) 


where O(dz*) denotes terms proportional to fourth power of dx. Similar terms 
with daz /dx2 and da3/6x3 will be given by contributions of the other faces so 
that for the whole parallelepiped whose volume dV = dz,dxz2dzx3 we have 
Equation: 


If we let the volume shrink to zero we have 


Equation: 
jm aq | ede Vea 


If a is a flux, then the surface integral is the net flux of a out of the volume. In 
particular, let a be the fluid velocity, which can be thought as a volumetric flux. 
Then the divergence of velocity is the volumetric expansion per unit volume. A 
vector field with identically zero divergence is called solenoidal. An 
incompressible fluid has a solenoidal velocity field. If the flux field of a certain 
property is solenoidal there is no generation of that property within the field, 
for all that flows into an infinitesimal element flows out again. 


If a is the gradient of a scalar function V@ , its divergence is called the 
Laplacian of ¢. 
Equation: 


7d Oo O7¢ 
V6 = Ve(Vd) = 4 = — + — + — 
? RD ee Oxi Ong One 


A function that satisfies Laplace's equation V7¢ = 0 is called a potential 
function or a harmonic function. An irrotational, incompressible flow field has 
a velocity that is the gradient of a flow potential. Also, the steady-state 
temperature field in a homogeneous solid and the steady state pressure 
distribution of a single fluid phase flowing in porous media are solutions of 
Laplace's equation. In two dimensions the solutions of Laplace's equation can 
be found through the use of complex variables. 


If A is a tensor, the notation divA or V e A is sometimes used for the vector 
Aj;;,;. The index notation is preferred for tensors. 


The curl of a vector field 


The symbolic vector product or cross product of the vector operator V anda 
vector field a(x) is called the curl of the vector field. It is the vector 
Equation: 


V xX a=curla = €j;% Ak; €(;) 


That this definition is a combination of the previously definitions for the V 
operator and the cross product can be seen be carrying out the operations. 


Equation: 


Vxa = (ew er + &2) ey FG) a) x (e(1) a1 + €(@) a2 + e(g) as) 
_ Oa 6) 0 Oa Oa O 
= ($22 — 5) e(1) + ($22 - se ) e(2) + ($2 = 5 ) e(3) 
also 
0a; 
VX a = Eijk Gy- C(k) = ijk @j,i C(k) = Ekji Vj,k C(’) = Ejki h,j C(;) 


= Cijk Qk,j e()Q. E. D. 


The connection between the curl of a vector field and the rotation of the vector 
field (it is called rot a in some older texts) can be illustrated by calculating the 
circulation of the vector field around a closed curve. Consider a elementary 
rectangle in the 023 plane normal to 01 with one comer P at (x1, %2, x3) and 
the diagonally opposite one @ at (11, 2 + dx2, x3 + dx3) as shown in Fig. 
3.8. We wish to calculate the line integral or circulation around this elementary 
closed curve of a e tds, where t is the unit tangent to the curve. Now the line 
through P parallel to 03 has tangent —e,3) and the parallel side through Q has 
tangent e(3), and each is of length dx3. Accordingly, they contribute to a e tds 


an amount 
Equation: 


O 
[as (x14, Lo+ dx, £3) — a3 (x1, £2, x3)| dx3 = = dx4 dx3 + O (dx*) 
2 


Similarly, from the other two sides, there is a contribution, 


Equation: 


Oa 
os dx3dx2 +O (dz*) 


Thus writing dA = dx2 dx3, we have 


Equation: 
1 0az3 Oa 
—— tds = | —— — —— O(d 
dA fas : & set ) + ide) 


023 
and in the limit 
Equation: 
1 Oaz Oa» 
lim —— tds = | -—- —- -— ] =(V aad oi 
vee dA fas ; (2 a | ( Re ( eae 
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The suffix 023 has been put on the integral sign to show that the line integral in 
a 023 plane, and the last equation shows that the circulation in the O23 plane is 
equal to the component of the curl in the O1 direction. This correspondence 
between the curl and circulation gives physical meaning the curl of a vector 
field. It is a measure of the circulation or rotation of the motion. There is a 
direction associated with circulation, rotation, and curl. If the circulation 
around a closed curve is in the direction of the closed fingers of the right hand, 
then the curl is in the direction of the thumb. 


A vector field a for which V x a = 0 is called irrotational because the 
circulation about any closed curve vanishes. 
Equation: 


Ve(Vxa)=0 

V x (Vd) =0 

V x (ga) = (Vb) xa+GVxa 

V x (ax b) =a(Veb)— b(Vea)+ (be V)a— (ae V)b 
Vx(Vxa)=V(Vea)—V2a 


The first of these states that if a vector field b is the curl of a vector, i.e., 

b = V x a then the vector field b is solenoidal. The second states that if a 
vector field b is equal to the gradient of a scalar, i.e., b = Vo then the vector 
field b is irrotational. The last identity has the Laplacian operator V2? = Ve V 
operating on a vector. The result is a vector whose components are equal to the 
Laplacian of the components, if the coordinates are Cartesian. This may not be 
the case in curvilinear coordinates. 


Green's theorem and some of its variants 


The divergence theorem, also called the Gauss' theorem, or Green's theorem 
equates the volume integral of the divergence of a vector field a(x) to the 
surface integral of the normal component of the vector. 


Equation: 
ia Veadv = fhaends 
V Ss 


During the discussion of the divergence of a vector field, we showed that the 
above relation holds for an infinitesimal volume. Suppose now that a 
macroscopic volume of space is a composite of the infinitesimal regions. The 
total volume integral is the sum of the infinitesimal volume integrals. However, 
the contribution of a e ndS' from the touching faces of two adjacent elements 
of volume are equal in magnitude but opposite in sign since the outward 
normal points in opposite directions. Thus in a summation of a e ndS, the only 
terms that survive are those on the outer surface S, i.e., the surface integral is 
over the exterior surfaces of the macroscopic region. Q.E.D. 


If a= Vo we have 
Equation: 


I V’ddV = II VeVddV = pris Pr 


where 6¢/6n denotes the derivative in the direction of the outward normal. If 
the scalar @ is temperature, this equation says that at steady-state, the integral 
of the net sources of heat in the volume is equal to the flux across the external 
surfaces. 


Stokes' theorem 


On a surface let C’ be the curve or a finite number of curves forming the 
complete boundary of an area S. We assume that the surface is two-sided and 
that S can be resolved into a finite number of regular elements. Choose a 
positive side of S and let the positive direction along C’ be that in which an 
observer on the positive side must move along the boundary if he is to have the 
area S always on his left. At each regular point on the surface let n be the unit 
normal drawn toward the positive side. Let a and its first derivatives be 
continuous on S. Stokes' theorem states that the circulation around a closed 
curve is equal to the surface integral of the normal component of the curl. 


Equation: 
fastds = [for a)endS 
C Ss 


We showed earlier the circulation around an infinitesimal, closed curve was 
equal to the normal component of the curl multiplied by the area of the 
enclosed surface. We will extend the earlier result for an infinitesimal closed 
curve enclosing an infinitesimal surface to a macroscopic curve and surface. 
The macroscopic surface will be subdivided into a composite of many 
infinitesimal regions where the earlier result apply. The summation of the 
normal component of the curl multiplied by the area of the element is equal to 


the surface integral of the normal component of the curl. However, the quantity 
a e tds from the touching sides of two adjacent surface elements have equal 
magnitude but opposite sign since the direction of the line integrals are in 
opposite directions. Thus in the summation of the circulations, the only terms 
that survive are the contribution of the external bounding curve, i.e., the 
circulation is around the exterior curve C bounding the surface S. Q.E.D. 


The classification and representation of vector fields 


We mentioned earlier that a solenoidal vector field is one where V e a = 0 
everywhere and an irrotational vector field is one where V x a = 0 
everywhere. A vector field that is the gradient of a scalar a = Vo is 
irrotational. If a vector field is both irrotational and solenoidal it is the gradient 
of a harmonic function, where V e (Vd) = V?¢ = 0. It can be proven that if 
a vector field is both irrotational and solenoidal, it is uniquely determined in a 
volume V if it is specified over S, the surface of V. There other types of 
named vector fields are discussed by Aris. 


Irrotational vector fields 


The vector field a is irrotational if its curl vanishes everywhere. By Stokes' 
theorem the circulation around any closed curve also vanishes. Also, an 
irrotational vector field can be expressed as the gradient of a scalar. 
Equation: 


Vxa-—0 


$ aetds=0 4 an irrotational field 


a= Vo 


The velocity field of motions where the viscous effects are insignificant 
compared to inertial effects and the flow is initially irrotational can be 
approximated as an irrotational velocity field. 


Solenoidal vector fields 


A solenoidal vector field is defined as one in which the divergence vanishes. 
This implies that the flux across a closed surface must also vanish. A vector 
identity states that the divergence of the curl of a vector is zero. Thus a 
continuously differentiable solenoidal vector field has the following three 
equivalent characteristics. 

Equation: 


Vea-—0 


Che endS=0 4 a solenoidal field 


a=VxA 


The velocity field of motions where the effects of compressibility are 
insignificant can be approximated as a solenoidal vector field. The surface 
integral of velocity vanishing over any closed surface means that the net 
volumetric flow across closed surfaces is zero. Incompressible flow fields can 
be expressed as the curl of a vector potential. Two-dimensional, incompressible 
flows have only one nonzero component of the vector potential and this is 
identified as the stream function. 


Helmholtz’ representation 


We found that an irrotational vector is the gradient of a scalar potential and a 
solenoidal vector is the curl of a vector potential. Here we show that any vector 
field with sufficient continuity is divisible into irrotational and solenoidal parts, 
and so is expressible in terms of a scalar and a vector potential. The 
fundamental problem in the analysis of a vector field is the determination of 
these potentials and their expression in terms of the essential characteristics of 
the vector, namely divergence, curl, discontinuities, and boundary values. For 
when the potentials are known the vector itself can be determined by 
differentiation. The following analysis is taken from H. B. Phillips, Vector 
Analysis, John Wiley & Sons, 1933. The following nomenclature will differ 
somewhat in that the vector is expressed as the negative of the gradient of a 
scalar. Also, the vector field of interest will be denoted as mathbf F’.. Bold face 


capital letters will also be used for other vector quantities. Also the equations 
have the 47 factor of electromagnetism in mks units rather than the factors €, 
and €, c? of the SI units. 


Let V be a region of space where the vector field F’ has piecewise continuous 
second derivatives, 5; be surfaces of discontinuity of F’, and S' be the 
bounding surface of V. The Helmholtz's theorem states that F’ can be 
expressed in terms of the potentials. 

Equation: 


F=-V6+VxA 


where 


bmp) = J{PEe+ [fees — gape 
V Sy Ss 
Aiee\= Iv , ppddS _ icp nxFas 
(ep) = JI J ee 


VeF=—-V*¢=4np 

Vx F=V~x (Vx A) =4aI 

ne AF = 470 

n x AF = 47J 

r = |xp — xg|wherexgiscoordinateofintegrand 


The vectors I and J are not arbitrary. They are subject to the equation of 
continuity V e A = 0. The effect of this condition is to make I and J behave 
like space and surface currents of something which is nowhere created or 
destroyed. In the electromagnetic field they usually represent currents of 
electricity. In hydrodynamic fields they represent vorticity. If A is everywhere 
solenoidal, the following three equations must then be everywhere satisfied. 
Equation: 


VelIl=0 


neAIl+VeJ=0 
nx AJ etds = 0 


Considering I and J as representing currents, the first equation expresses that 
the amount which flows out of a small region is equal to the amount which 
flows in. The second equation expresses that the total flow from a portion of a 
conducting surface into space and along the surface is zero. The third equation 
expresses that the flow from one sub-region across a curve on a conducting 
surface is equal to the flow into the adjacent sub-region. These three equations 
thus express that I and J, considered as space and surface currents, represent a 
flow of something which is conserved. For I and J to have this property the 
above discussion shows it is necessary and sufficient that A be everywhere 
solenoidal. In hydrodynamics, the field I corresponds to vorticity and it clearly 
is solenoidal because it is the curl of velocity. 


Vector and scalar potential 


The previous section showed that a vector field can be determined from the 
divergence and curl of the vector field and the values on surfaces of 
discontinuities and bounding surfaces. The integral equations are useful for 
developing analytical solutions for simple systems. However, in 
hydrodynamics the vorticity is generally an unknown quantity. Thus it is useful 
to express the potentials as differential equations that are solved simultaneously 
for the potentials and vorticity. The differential equations are derived by 
substituting the potentials into the expressions for the divergence and curl. 
Equation: 


V*6 = —4rp 
V2A = —4rI 


In two-dimensional vector fields the vector potential A and the vector I has a 
nonzero component only in the third direction. In hydrodynamics the nonzero 
component of the vector potential is the stream function 47I corresponds to the 
vorticity, which has only one nonzero component. 


Assignment 3.1: Particle velocity and acceleration 
Suppose a particle fixed on the surface of a steady-rotating sphere with radius, 
R, has a constant speed, |v| = v. 


a. Show that its acceleration is perpendicular to its velocity. 

b. Show that the acceleration has a radial (from the center of the sphere) 
component a, = —v?/R. 

c. Let the magnitude of the angular velocity be w. Express the centrifugal 
acceleration (direction perpendicular to the axis of rotation), a in terms of 
R, w and the angle of particle from the axis of rotation. 


Assignment 3.2: Differential area and volume 


a. Express the differential of area in term of the differential of the surface 
coordinates for a spherical surface using the spherical polar coordinate 
system. 

b. Obtain the differential volume elements in cylindrical and spherical polars 
by the Jacobian and check with a simple geometrical picture. 


Assignment 3.3: Differential operators 


a. Derive the expression for the Laplacian of a scalar in Cartesian 
coordinates from the definition of the gradient and divergence. 

b. Prove that: Ve (da) = Vb ea+ o(Vea) 

c. Let x be the Cartesian coordinates of points in space and r = |x| . 
Calculate the divergence and curl of x and the gradient and Laplacian of r 
and 1/r. Note any singularities. 

d. Prove the identities involving the curl operator. 

e. Suppose a rigid body has the velocity field v = vq + w x (x — X,). 
Show that the curl of this velocity field is V x v = 2w. 


The Kinematics of Fluid Motion 
Topics Covered in this Chapter 


e Particle paths and material derivatives 

e Streamlines 

e Streaklines 

e Dilatation 

e Reynolds' transport theorem 

¢ Conservation of mass and the equation of continuity 
¢ Deformation and rate of strain 

e Physical interpretation of the deformation tensor 

e Principal axis of deformation 

e Vorticity, vortex lines, and tubes 


Reading assignment: Chapter 4 of Aris 


Kinematics is the study of motion without regard to the forces that bring about the motion. 
Already, we have described how rigid body motion is described by its translation and rotation. 
Also, the divergence and curl of the field and values on boundaries can describe a vector field. 
Here we will consider the motion of a fluid as microscopic or macroscopic bodies that translate, 
rotate, and deform with time. We treat fluids as a continuum such that the fluid identified to be at 
a specific point in space at one time with neighboring fluid will be at another specific point in 
space at a later time with the same neighbors, with the exception of certain bifurcations. This 
identification of the fluid occupying a point in space requires that the motion is deterministic 
rather than stochastic, i.e., random motions such as diffusion and turbulence are not described. 
Central to the kinematics of fluid motion is the concept of convection or following the motion of 
a "particle" of fluid. 


Particle paths and material derivatives 


Fluid motion will be described as the motion of a "particle" that occupies a point in space. At 
some time, say t = 0, a fluid particle is at a position € = (£1, 2, €3) and at a later time the same 
particle is at a position x. The motion of the particle that occupied this original position is 
described as follows. 

Equation: 


x= (6,1) 08 wy —"e (61;60}£3; 8) 


The initial coordinates € of a particle will be referred to as the material coordinates of the 
particles and, when convenient, the particle itself may be called the particle € . The terms 
convected and Lagrangian coordinates are also used. The spatial coordinates x of the particle 
may be referred to as its position or place. It will be assumed that the motion is continuous, single 
valued and the previous equation can be inverted to give the initial position or material 
coordinates of the particle which is at any position x at time ¢; i.e., 

Equation: 


é= E (x, t) or bi _ bi (x1, £2, £3, t) 


are also continuous and single valued. Physically this means that a continuous arc of particles 
does not break up during the motion or that the particles in the neighborhood of a given particle 
continue in its neighborhood during the motion. The single valuedness of the equations mean that 
a particle cannot split up and occupy two places nor can two distinct particles occupy the same 
place. Exceptions to these requirements may be allowed on a finite number of singular surfaces, 
lines or points, as for example a fluid divides around an obstacle. It is shown in Appendix B that 
a necessary and sufficient condition for the inverse functions to exist is that the Jacobian 
Equation: 


O(x1, £2, £3) 


<= O(E1, €2, €3) 


should not vanish. 


The transformation x = x(, t) may be looked at as the parametric equation of a curve in space 
with ¢ as the parameter. The curve goes through the point € , corresponding to the parameter 

t = 0, and these curves are the particle paths. Any property of the fluid may be followed along 
the particle path. For example, we may be given the density in the neighborhood of a particle as a 
function p(€, t), meaning that for any prescribed particle € we have the density as a function of 
time, that is, the density that an observer riding on the particle would see. (Position itself is a 
"property" in this general sense so that the equations of the particle path are of this form.) This 
material description of the change of some property, say 3(€, t), can be changed to a spatial 
description 3(x, t). 

Equation: 


3(x, t) = 3[E(x, t), ¢] 


Physically this says that the value of the property at position x at time ¢ is the value appropriate 
to the particle that is at x at time ¢. Conversely, the material description can be derived from the 
spatial one. 
Equation: 


S(g,t) = F[E(x, t), ¢] 
meaning that the value as seen by the particle at time t is the value at the position it occupies at 


that time. 
Equation: 


O 6) ere : : ; 
aE i= derivative with respect to time keeping x constant 
x 


and 
Equation: 


D O tae . ; ; 
— =| —] = derivative with respect to time keeping € constant 
Dt Ot é 


Thus 03 /0t is the rate of change of 3 as observed at a fixed point x, whereas D7/ Dt is the rate 
of changed as observed when moving with the particle, i.e., for a fixed value of € . The latter we 
call the material derivative. It is also called the convected, convective, or substantial derivative 
and often denoted by D/ Dt. In particular the material derivative of the position of a particle is it 
velocity. Thus putting J = z;, we have 


Equation: 
oa = pH, Yt »§2, €3,¢ 
De ee (£1, €2, €3, t) 
or 
Equation: 
_ Dx 
Dt 


This allows us to establish a connection between the two derivatives, for 
Equation: 


D3 0. mec 
Dp = gi Gt= 3 Ie Got 


_ (8) , 8 (dx; 
A 


a, oO 


Streamlines 


We now have a formal definition of the velocity field as a material derivative of the position of a 
particle. 


Equation: 


ens Dr) 2 ont 


The field line lines of the velocity field are called streamlines; they are the solutions of the three 
simultaneous equations 


Equation: 


where s is a parameter along the streamline. This parameter s is not to be confused with the time, 
for in the above equation t is held fixed while the equations are integrated, and the resulting 
curves are the streamlines at the instantt. These may vary from instant to instant and in general 
will not coincide with the particle paths. 


To obtain the particle paths from the velocity field we have to follow the motion of each particle. 
This means that we have to solve the differential equations 
Equation: 


subject to x = € at t = 0. Time is the parameter along the particle path. Thus the particle path is 
the trajectory taken by a particle. 


The flow is called steady if the velocity components are independent of time. For steady flows, 
the parameter s along the streamlines may be taken to be ¢ and the streamlines and particle paths 
will coincide. The converse does not follow as there are unsteady flows for which the streamlines 
and particle paths coincide. 


If C is a closed curve in the region of flow, the streamlines through every point of C’ generate a 
surface known as a stream tube. Let S be a surface with C’ as the bounding curve, then 


Equation: 
/ vendS 


is known as the strength of the stream tube at its cross-section S. 


The acceleration or the rate of change of velocity is defined as 


Equation: 
Dv Ov 
— —_—_— V 
a Di At +(veV)v 
Ov; he Ov; 
Qa = VU; 
ot : Ox j 


Notice that in steady flow this does not vanish but reduces to 
Equation: 


a= (ve V)v for steady flow. 


Even in steady flow other than a constant translation, a fluid particle will accelerate if it changes 
direction to go around an obstacle or if it increases its speed to pass through a constriction. 


Streaklines 


The name streakline is applied to the curve traced out by a plume of smoke or dye, which is 
continuously injected at a fixed point but does not diffuse. Thus at time ¢ the streakline through a 
fixed-point y is a curve going from y to x(y, t), the position reached by the particle which was at 
y at time t = 0. A particle is on the streakline if it passed the fixed-point y at some time between 
0 and ¢. If this time was t’, then the material coordinates of the particle would be given by 

£ = E(y, t'). However, at time t this particle is at x = x(€, t) so that the equation of the 
streakline at time ¢ is given by 

Equation: 


“= xE(y,t’), €], 


where the parameter ¢’ along it lies in the interval 0 < t’ < t. If we regard the motion as having 
been proceeding for all time, then the origin of time is arbitrary and t' can take negative values 
—oo <dt/ <t. 


The flow field illustrated in 4.13 by Aris is assigned as an exercise. 


Dilatation 


We noticed earlier that if the coordinate system is changed from coordinates € to coordinates z, 
then the element of volume changes by the formula 
Equation: 


O(x1, £2, z3) 


dV = 
(1, £9, £3) 


d&, d&, d€3 = J dVo 


If we think of € as the material coordinates, they are the Cartesian coordinates at t = 0, so that 
d€,d&2dé&3 is the volume dV, of an elementary rectangular parallelepiped. Consider this 
elementary parallelepiped about a given point € at the initial instant. By the motion this 
parallelepiped is moved and distorted but because the motion is continuous it cannot break up and 
so at some time t is some neighborhood of the point x = x(€,t). By the above equation, its 
volume is dV = J dV, and hence 

Equation: 


V . . $5 aS at 
5 = ratio of an elementary material volume to its initial volume. 


It is called the dilation or expansion. The assumption that x = x(€, t) can be inverted to give 
€ = €(x, t), and vice versa, is equivalent to requiring that neither J nor J~ 1 vanish. Thus, 
Equation: 


0<J<o 


We can now ask how the dilation changes as we follow the motion. To answer this we calculate 
the material derivative DJ / Dt. However, 


Equation: 
Oxy Ot, Ot, 
Of: Of 0&3 
Fetes Ox; Ou; OxR Or, 0x2 Oxy 
OE, Ok OE, | OH A OE 
Ox3 0x3 0x3 
Of: Of 0&3 
Now 
Equation: 


0g; og; Dt = 0; 


for D/ Dt is differentiation with € constant so that the order can be interchanged. Now if we 
regard v; as a function of x1, £2, £3, 
Equation: 

Ov; Ov; Om 

Of; 08; 08; 


The above relation can now be applied to differentiation of the Jacobian. 


Equation: 
DJ D [ 0x; \ 0X; Oxz, Ox; D (02; \ Oxy Ox; Or; D ( Oxy 
De > ee (se) OE, Oe; |“ "* OE, Dt (se) G6. OE Ob; Di (se) 
sg OUI Os ONO yes, OE UN Orn Ot i, a ORs UP OUn Oto, 
"Om O& O&s Es OE, O%m Ob, OEs |" OE, OE, Oxm OE; 
Ov, Ox; OX; Oxy, Ox; Ov; OX; Ox, uct Ox; Ov3 


rey + E45 +€ 
* On, Of; Of, O€3 |” OE, Oxy Oly OE, |” OE, Oka Azz 
Ovi Ov2 O0v3 
= J+ J J 
Ox Ox2 uy 023 


= (Vev)J 


where we made use of the property of the determinant that the determinant of a matrix with 
repeated rows is zero. Thus, 
Equation: 


D(in J) 


=Vev 
Dt 


We thus have an important physical meaning for the divergence of the velocity field. It is the 
relative rate of dilation following a particle path. It is evident that for an incompressible fluid 
motion, 

Equation: 


V ev =0 for incompressible fluid motion. 


Use of a stream function to satisfy the mass-conservation equation (Batchelor, 
1967) 


In the cases of flow of an incompressible fluid, and of steady flow of a compressible fluid, the 
mass-conservation equation reduces to the statement that a vector divergence is zero, the 
divergences being of u and pu respectively. If we impose the further restriction that the flow 
field either is two-dimensional or has axial symmetry, this vector divergence is the sum of only 
two derivatives, and the mass-conservation equation can then be regarded as defining a scalar 
function from which the components of u or pu are obtained by differentiation. The procedure 
will be described here for the case of an incompressible fluid. 


Assume first that the motion is two-dimensional, so that u = (u, v, 0) and u and v are 
independent of z. The mass-conservation equation for an incompressible fluid then has the form 
Equation: 


du, av _, 
Ox Oy 


from which it follows that wy — vdz is an exact differential, equal to dw say. Then 
Equation: 


and the unknown scalar function (2, y, t) is defined by 
Equation: 


U—tbo= | (udy—v dr), 


where w, is a constant and the line integral is taken along an arbitrary curve joining some 
reference point O to the point P with co-ordinates z, y. In this way we have used the mass- 
conservation equation to replace the two dependent variables u, v by the single dependent 
variable ~ , which is a very valuable simplification in many cases of two-dimensional flow. 


Calculation of the flux of fluid volume across a curve 
joining the reference point O to the point P with co- 
ordinates (x, y) 


The physical content of the above argument also proves to be of interest. The flux of fluid volume 
across a curve joining the points 0 and P in the (a, y)-plane (by which is meant the flux across 
the open surface swept out by translating this curve through unit distance in the z-direction), the 
flux being reckoned positive when it is in the anti-clockwise sense about P, is given exactly by 
the right-hand side of (2.2.8) (see figure 2.2.1). Now the flux of volume across the closed curve 
formed from any two different paths joining 0 to P is necessarily zero when the region between 
the two paths is wholly occupied by incompressible flow. The flux represented by the integral in 
(2.2.8) is therefore independent of the choice of the path joining O to P, provided it is one of a 
set of paths of which any two enclose only incompressible fluid, and therefore defines a function 
of the position of P, which we have written as w — wo. 


Since the flux of volume across any curve joining two points is equal to the difference between 
the values of ~ at these two points, it follows that w= is constant along a streamline, as is also 
apparent from (2.2.7) and the equation (2.1.1) that defines the streamlines. ~ is termed the stream 
function, and is associated (in this case of two-dimensional flow) with the name of Lagrange. The 
function can also be regarded as the only non-zero component of a 'vector potential’ for u 
(analogous to the vector potential of the magnetic induction, which also is a solenoidal vector, in 
electromagnetic field theory), since (2.2.7) can be written as 

Equation: 


u=VxA, A=(0,0,7) 


It is common practice in fluid mechanics to provide a picture of a flow field by drawing various 
streamlines, and if these lines are chosen so that the two values of ~ on every pair of neighboring 
streamlines differ by the same amount, € say, the eye is able to perceive the way in which the 
velocity magnitude q, as well as its direction, varies over the field, since 

Equation: 


q = €//(distance between neighboring streamlines). 


Examples of families of streamlines describing two-dimensional flow fields, with equal intervals 
in = between all pairs of neighboring streamlines, will be found in figures 2.6.2 and 2.7.2. 


Expressions for the velocity components parallel to any orthogonal coordinate lines in terms of w 
may be obtained readily, either by the use of (2.2.9) or with the aid of the relation between w and 
the volume flux 'between two points’. For flow referred to polar co-ordinates (7,9), we find, by 
evaluating the flux between pairs of neighboring points on the r- and 6-co-ordinate lines and 
equating it to the corresponding increments in w (allowance being made for the signs in the 
manner required by (2.2.8)), that 

Equation: 


The reader may find useful the general rule for two-dimensional flow, that differentiation of ~ in 
a certain direction gives the velocity component 90 degree in the clockwise sense from that 
direction. 


Finally, for this case of two-dimensional flow of incompressible fluid, we should note the 
possibility that 7 is a many-valued function of position. For suppose that across some closed 
inner geometrical boundary there is a net volume flux m; this flux might be due to an effective 
creation of fluid within the inner boundary (as when a tube discharges fluid into this region) or to 
change of volume of the part of the enclosed region not occupied by the fluid (as when a gaseous 
cavity surrounded by water expands or contracts). 


If now we choose two different paths joining the two points 0 and P which together make up a 
closed curve enclosing the inner boundary, the fluxes of volume across the two joining curves 
differ by an amount m (or, more generally, by pm, where p is the number of times the combined 
closed curve passes round the inner boundary). The value of 7 — a, at the point P thus depends 
on the choice of path joining it to the reference point 0, and may take any one of a number of 
values differing by multiples of m. This kind of many-valuedness of a scalar function related to 
the velocity distribution in a region which is not singly-connected will be described more fully in 
§ 2.8. It is not confined to two-dimensional flow, although that is the context in which it occurs 
most often. 


If now the flow has symmetry about an axis, the mass-conservation equation for an 
incompressible fluid takes the form 
Equation: 


_ Ou x 1 O(ov) 


=0 
Ox o Oo 


in terms of cylindrical co-ordinates (x, 0, y) with corresponding velocity components (u, v, w), 
the axis of symmetry being the line o = 0 . This relation ensures that cwdo — ovdz is an exact 
differential, equal to dy say. Then 

Equation: 


and the function (2, a, t) is defined by 
Equation: 


y-%o= [o(uae - vdz), 


where the line integral is taken along an arbitrary curve in an axial plane joining some reference 
point 0 to the point P with co-ordinates (x, o). It will be noticed that the azimuthal component of 
velocity w does not enter into the mass-conservation equation in a flow field with axial symmetry 
and cannot be obtained from ~w. 


Again it is possible to interpret 7 both as a measure of volume flux and as one component of a 
vector potential. The flux of fluid volume across the surface formed by rotating an arbitrary curve 
joining 0 to P in an axial plane, about the axis of symmetry, the flux again being reckoned as 
positive when it is in the anti-clockwise sense about P, is 27-times the right-hand side of 
(2.2.12). Lines in an axial plane on which w is constant are everywhere parallel to the vector (uw, v 
, 0), and can be described as “streamlines of the flow in an axial plane’. 7 is here termed the 
Stokes stream function. A sketch of lines on which w is constant, with the same increment in w 
between all pairs of neighboring lines (see figure 2.5.2 for an example), does not give quite as 
direct an impression of the distribution of velocity magnitude here as in two-dimensional flow, 
owing to the occurrence of the factor 1/o in the expressions for wu and v in (2.2.11). The relations 
(2.2. 11) are readily seen to be equivalent to 

Equation: 


u=VxA, A,=v/o 


the components of the vector potential A referred to cylindrical co-ordinate lines here being 
independent of the azimuthal angle yp. 


The relation between and the volume flux 'between two points' may be used to obtain expressions 
for the velocity components referred to other orthogonal systems of co-ordinates in terms of w. 
For instance, for flow with axial symmetry referred to spherical polar co-ordinates (r, 0, y), we 
find, by evaluating the flux between pairs of neighboring points on the r- and 0 -coordinate lines 


and equating it to 27 times the corresponding increments in w (allowance being made for the 
signs in the manner required by (2.2.12)), that 
Equation: 
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With this co-ordinate system, the vector potential for the velocity has the azimuthal component 
Equation: 


Exercise 


At time to the position of a material element of fluid has Cartesian coordinates (a, b, c) and the 
density of the fluid is p,. At a subsequent time t the position coordinates and density of the 
element are (X, Y, Z) and p . Show that with this Lagrangian specification of the flow field the 
equation of mass conservation is 

Equation: 


Reynolds’ transport theorem 


An important kinematical theorem can be derived from the expression for the material derivative 
of the Jacobian. It is due to Reynolds and concerns the rate of change not of an infinitesimal 
element of volume but any volume integral. Let 3(x, t) be any function and V(t) be a closed 
volume moving with the fluid, that is consisting of the same fluid particles. Then 


Equation: 
F(t) = [[ [reoav 


v(t) 


is a function of ¢ that can be calculated. We are interested in its material derivative DF'/ Dt. Now 
the integral is over the varying volume V(t) so we cannot take the differentiation through the 
integral sign. If, however, the integration were with respect to a volume in €-space it would be 
possible to interchange differentiation and integration since D/Dt is differentiation with respect 
to t keeping € constant. The transformation x = x(,t) ,dV = JdV, allows us to do just this, 
for V(t) has been defined as a moving material volume and so come from some fixed volume V, 
at t = 0. Thus 


Equation: 


& [from — § |ffomsore 
- [Ge 2) 
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v(t) 


Since D/Dt = (0/0t) + v e V we can express this formula into a number of different forms. 
Substituting for the material derivative and collecting the gradient terms gives 


Equation: 
< [[ [resoav = ie bveva43(Vew) av 


v(t) v(t) 
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v(t) 


Now applying Green's theorem to the second integral we have 


Equation: 
off fresav= [fom av + $ ov) enas 


v(t) S(t) 


where S(t) is the bounding surface of V(t). This admits of an immediate physical picture for it 
says that the rate of change of the integral of 3 within the moving volume is the integral of the 
rate of change of 3 at a point plus the net flow of over the bounding surface. 3 can be any scalar 
or tensor component, so that this is a kinematical result of wide application. It is going to be the 
basis for the conservation of mass, momentum, energy, and species. This approach to the 
conservation equations differs from the approach taken by Bird, Stewart, and Lightfoot. They 
perform a balance on a fixed volume of space and explicitly account for the convective flux 
across the boundaries. 


Conservation of mass and the equation of continuity 


Although the idea of mass is not a kinematical one, it is convenient to introduce it here and to 
obtain the continuity equation. Let p(x, t) be the mass per unit volume of a homogeneous fluid at 
position x and time t. Then the mass of any finite material volume V(t) is 


Equation: 
n= [[ focsav. 
v(t) 


If V is a material volume, that is, if it is composed of the same particles, and there are no sources 
or sinks in the medium we take it as a principle that the mass does not change. By inserting J = p 
in Reynolds' transport theorem we have 


Equation: 
dm Dp 
a6 [[[{Browenhay 
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This is true for an arbitrary material volume and hence the integrand itself must vanish 
everywhere. It follows that 
Equation: 


which is the equation of continuity. 


A fluid for which the density p is constant is called incompressible. In this case the equation of 
continuity becomes 
Equation: 


Vev =O incompressible flow 


and the motion is isochoric or the velocity field solenoidal. 


Combining the equation of continuity with Reynolds’ transport theorem for a function J = pF 
we have 
Equation: 


i [fore = [ff ion <orceeo)w 


u(t) v(t) 


Liha ees)? 
|ffoBea 


v(t) 
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This equation is useful for deriving the conservation equation of a quantity that is expressed as 
specific to a unit of mass, e.g., specific internal energy and species mass fraction. 


Deformation and rate of strain 


The motion of fluids differs from that of rigid bodies in the deformation or strain that occurs with 
motion. Material coordinates give a reference frame for this deformation or strain. 


Consider two nearby points P and Q with material coordinates € and € + dé. At time t they are to 
be found at x(€, t) and x(€ + dé, t). Now 
Equation: 


Ox; 
fe ge 40d") 
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where O (d?) represents terms of order dé” and higher which will be neglected from this point 
onward. Thus the small displacement vector dé has now become 


Equation: 
= x (€ + dé, t) = x (€, t) 
where 
Equation: 
Ox; 
dx; = d€;. 
56, 


It is clear from the quotient rule (since d€ is arbitrary) that the nine quantities Ox; /O€,; are the 
components of a tensor. It may be called the displacement gradient tensor and is basic to the 
theories of elasticity and rheology. For fluid motion, its material derivative is of more direct 
application and we will concentrate on this. 


If v = Dx/Dt is the velocity, the relative velocity of two particles € and € + dé has components 
Equation: 


However, by inverting the above relation, we have 
Equation: 

_ Ov; OL _— Ov; 
(Oy Ox; 7? On; 


dv; 


expressing the relative velocity in terms of current position. Again it is evident that the 
(Ov; /Ox is) are components of a tensor, the velocity gradient tensor, for which we need to obtain a 
sound physical feeling. 


We first observe that if the motion is a rigid body translation with a constant velocity u, 
Equation: 


x=€+ut 


and the velocity gradient tensor vanishes identically. Secondly, the velocity gradient tensor can be 
written as the sum of symmetric and antisymmetric parts, 


Equation: 
Ovi (ge +5) +3 ($s — 
OF | DOG: * OD; 2\ Ox; Ox; 
= ey t 2; 
or 
Vv = e+ 


We have seen that a relative velocity dv; related to the relative position dx; by an antisymmetric 
tensor 2;;, i.e., dv; = 2;; dx;, represents a rigid body rotation with angular velocity 
w = —vec®1 . In this case 


Equation: 
1 1 Ov 
Wi = ~ 5 Siskin = sti a 
or 
1 
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Thus the antisymmetric part of the velocity gradient tensor corresponds to rigid body rotation, 
and, if the motion is a rigid one (composed of a translation plus a rotation), the symmetric part of 
the velocity gradient tensor will vanish. For this reason the tensor e;; is called the deformation or 


rate of strain tensor and its vanishing is necessary and sufficient for the motion to be without 
deformation, that is, rigid. 


Physical interpretation of the (rate of) deformation tensor 


The (rate of) deformation tensor is what distinguishes fluid motion from rigid body motion. 
Recall that a rigid body is one in which the relative distance between two points in the body does 
not change. We show here that the (rate of) deformation tensor describes the rate of change of the 
relative distance between two particles in a fluid. Also, it describes the rate of change of the angle 
between three particles in the fluid. 


First we will see how the distance between two material points change during the motion. The 
length of an infinitesimal line segment from P to Q is ds, where 
Equation: 


2 = . 2 — 
dx* = dx;dz; aE; Ob 
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now P and Q are the material particles € and € + d€ so that dé; and d€;, do not change during the 
motion. Also, recall that Dx; /Dt = v;. Thus 


Equation: 
D Ov; Ox; Ox; Ov; Ov; Ox; 
pees d 2 = 7 v7 7 v7 d d = y) v7 7 d d 
Be") = (See Be, + See Bey ) ele = 2p Fe te 
by symmetry. However, 
Equation: 
Ov; 7 Ou; Oxi 
ae, dé; = dvp= Be, dz; and BE, d&, = dz, 
since v = v|x(€)| and x = x(€) 
Thus 
Equation: 
5 Dy ) = (ds) 5, (ds) = 5 dz ;dx; = (e4; + N;;)dx dx; = e,jdzx dz; 
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by symmetry, or 
Equation: 


Now dz;/ds is the ith component of a unit vector in the direction of the segment PQ, so that this 
equation says that the rate of change of the length of the segment as a fraction of its length is 
related to its direction through the deformation tensor. 


In particular, if PQ is parallel to the coordinate axis 01 we have dx/ds = e€(1) and 
Equation: 
1 


ae Di) = e€1; in direction of 01 


Thus e 1, is the rate of longitudinal strain of an element parallel to the 01 axis. Similar 
interpretations apply to €g9 and €33. 


Now let's examine the angle between two line segments during the motion. Consider the segment, 
PQ and PR where PQ is the segment € + dé as before and PR is the segment € + dé’. If 6 is 
the angle between them and ds’ is the length of PR, we have from the scalar product, 

Equation: 


ds ds' cos 0 = dz;dz’, 


Differentiating with respect to time we have 
Equation: 


D(dz;dz') 

Dt 
— du;dz’, + dx;dv, 
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= 2e,;dz ;dz, 


since du; = (Ov; /Ox;)dzx’,. The i and j are dummy suffixes so we may interchange them in the 


first term on the right, then performing differentiation we have after dividing by dsds’ 
Equation: 


=z [dsds’ cos 8] = cos fares -- = ate} sin = 
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Now suppose that dz’ is parallel to the axis 01 and dz to the axis 02, so that 
(dz, /ds') and (dx; /ds) — 652 and O10 = m/2 . Then 
Equation: 


Thus €jg is to be interpreted as one-half the rate of decrease of the angle between two segments 
originally parallel to the 01 and 02 axes respectively. Similar interpretations are appropriate to e23 
and e€3). 


The fact that the rate of deformation tensor is linear in the velocity field has an important 
consequence. Since we may superimpose two velocity fields to form a third, it follows that the 
deformation tensor of this is the sum, of the deformation tensors of the fields from which it was 
superimposed. If v; = A (4)a;(no summation on 7) , we have a deformation which is the 
superposition of three stretching parallel to the three axis. However, if v1 = f (2), v2 = 0, 

v3 = 0 so that only nonzero component of the deformation tensor is e;2 = “f' (x2), the motion 
is one of pure shear in which elements parallel to the coordinate axis is not stretched at all. Note 
however that in pure stretching an element not parallel or perpendicular to the direction of 
stretching will suffer rotation. Likewise in pure shear an element not normal to or in the plane of 
shear will suffer stretching. 


Principal axis of deformation 


The rate of deformation tensor is a symmetric tensor and the principal axis of deformation can be 
found. They correspond to the eigenvalues of the matrix and the eigenvalues are the principal 
rates of strain. A set of particles that is originally on the surface of a sphere will be deformed to 
an ellipsoid whose axes are coincident with the principal axis. 


Vorticity, vortex lines, and tubes 


We have frequently reminders of rotating bodies of fluid such as tropical storms, hurricanes, 
tornadoes, dust devils, whirlpools, eddies in the flow behind objects, turbulence, and the vortex in 
draining bathtubs. The kinematics of these fluid motions is described by the vorticity. 


The antisymmetric part of the rate of strain tensor §2;; represents the local rotation, w;. Recall 
Equation: 


Wi = 5 aj Dj 


The curl of velocity is known as the vorticity, 
Equation: 


w=Vxv 


Thus the vorticity and the antisymmetric part of the rate of strain tensor is a measure of the 
rotation of the velocity field. An irrotational flow field is one in which the vorticity vanishes 
everywhere. The field lines of the vorticity field are called vortex lines and the surface generated 
by the vortex lines through a closed curve C is a vortex tube. The strength of a vortex tube is 
defined as the surface integral of the normal component. It is equal to the circulation around the 
closed curve C’ that bounds the cross-section S by Stokes' theorem. 


Equation: 
/ wendS = [fx v) ends 
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We observe that the strength of a vortex tube at any cross-section is the same, as w is a solenoidal 
vector. The surface integral of the normal component of a solenoidal vector vanishes over any 
closed surface. The surface integral on the surface of a vortex tube is zero because the sides are 
tangent to the vorticity vector field. Thus the surface integral across any cross-section must be 
equal in magnitude. The magnitude of the vorticity field can be visualized from the relative width 
of the vortex tubes in the same manner that the magnitude of the velocity field can be visualized 
by the width of the stream tubes. 


Because the strength of the tube does not vary with position along the tube, it follows that the 
vortex tubes are either closed, go to infinity or end on solid boundaries of rotating objects. In a 
real fluid satisfying the no-slip boundary condition, vortex lines must be tangential to the surface 
of a body at rest, except at isolated points of attachment and separation, because the normal 
component of vorticity vanishes on the stationary solid. 


When the vortex tube is immediately surrounded by irrotational fluid, it will be referred to as a 
vortex filament. A vortex filament is often just called a vortex, but we shall use this term to 
denote any finite volume of vorticity immersed in irrotational fluid. Of course, the vortex 
filament and the vortex require the fluid to be ideal (zero viscosity) to make strict sense, because 
viscosity diffuses vorticity, but they are useful approximations for real fluids of small viscosity. 


Helmholtz gave three laws of vortex motion in 1858. For the motion of an ideal (zero viscosity) 
barotropic (density is a single valued function of pressure) fluid under the action of conservative 
external body forces (gradient of a scalar), they can be expressed as follows: 


I. Fluid particles originally free of vorticity remain free of vorticity. 
II. Fluid particles on a vortex line at any instant will be on a vortex line at subsequent times. 
Alternatively, it can be said that vortex lines and tubes move with the fluid. 
II. The strength of a vortex tube does not vary with time during the motion of the fluid. 


The equations for the dynamics of vorticity will be developed later. 


P. G. Saffman, Vortex Dynamics, Cambridge University Press, 1992. C. Truesdell, The 
Kinematics of Vorticity, Indiana University Press, 1954. 


Assignment 4.1 

Plot the streamlines, particle paths, and streaklines of the flow field described in Sec. 4.13. Find 
the CHBE 501 web page. Download the files in CHBE501/Problems/lines. Execute lines with 
MATLAB. 


Assignment 4.2 

Execute the program deform in CHBE 501/Problems/deform and print the figures. It computes 
the particle paths for a patch of particles deforming in Couette flow, stagnation flow, and that of 
Sec 4.13. Look at the respective subroutine to determine the equations of the flow field. For each 
of these flow fields calculate: 


a. divergence 

b. curl 

c. rate of deformation tensor 

d. antisymmetric tensor 

e. Which fields are solenoidal or irrotational? 


Stress in Fluids 
Topics Covered in this Chapter 


e Cauchy's stress principle and the conservation of momentum 
e The stress tensor 

e The symmetry of the stress tensor 

e Hydrostatic pressure 

e Principal axes of stress and the notion of isotropy 

e The Stokesian fluid 

e Constitutive equations of the Stokesian fluid 

e The Newtonian fluid 

e Interpretation of the constants A and yu 


Reading assignment 
Chapter 1 in BSL 
Chapter 5 in Aris 


The only material property of the fluid we have so far discussed is the 
density. In the last chapter we introduced the rate of deformation or rate of 
strain tensor. The distinguishing characteristic between fluids and solids is 
that fluids can undergo unlimited deformation and yet maintain its integrity. 
The relation between the rate of deformation tensor and stress tensor is the 
mechanical constitutive equation of the material. An ideal fluid has a stress 
tensor that is independent of the rate of deformation, i.e., it has an isotropic 
component, which is identified as the pressure and has zero viscosity. 


Cauchy's stress principle and the conservation of momentum 


The forces acting on an element of a continuous medium may be of two 
kinds. External or body forces, such as gravitation or electromagnetic 
forces, can be regarded as reaching into the medium and acting throughout 
the volume. If the external force can be describes as the gradient of a scalar, 
the force is said to be conservative. Internal or contact forces are to be 
regarded as acting on an element of volume through its bounding surfaces. 
If an element of volume has an external-bounding surface, the forces there 
may be specified, e.g., when a constant pressure is applied over a free 
surface. If the element is internal, the resultant force is that exerted by the 


material outside the surface upon that inside. Let n be the unit outward 
normal at a point of the surface S and t,,,) the force per unit area exerted 
there by the material outside S. Then Cauchy's principle asserts that t(,,) is 
a function of the position x, the time ¢, and the orientation n of the surface 
element. Thus the total internal force exerted on the volume V through the 


bounding surface S is 
} / t(nydS. 


Equation: 
If f is the external force per unit mass (e.g. if 03 is vertical, gravitation will 
exert a force —ge,3) per unit mass or —pge,3) per unit volume, the total 


external force will be 
/ / / ofdV. 


Equation: 
UV 


The principle of the conservation of linear momentum asserts that the sum 
of these two forces equals the rate of change of linear momentum of the 
volume, i.e., 


Equation: 
D 


This is just a generalization of Newton's law of motion, which states that 
the rate of change of momentum of a particle is equal to the sum of forces 
acting on it. It has been extended to a volume that contains a number of 
particles. 


From the form of these integral relations we can deduce an important 
relation. Suppose V is a volume of a given shape with characteristic 
dimension d. Then the volume of V will be proportional to d° and the area 
of S to d’, with the proportionality constants depending only on the shape. 
Now let V shrink to a point but preserve its shape, then the volume 
integrals in the last equation will decrease as d® but the surface integral will 
decrease as d?. It follows that 


Equation: 
lim : t/,dS = 0 
0 (n) 


or, the stresses are locally in equilibrium. 


The stress tensor 


To elucidate the nature of the stress system at a point P we consider a small 
tetrahedron with three of its faces parallel to the coordinate planes through 
P and the fourth with normal n (see Fig. 5.1 of Aris). If dA is the area of 
the slanted face, the areas of the faces perpendicular to the coordinate axis 
Piis 

Equation: 


The outward normals to these faces are —e,;) and we may denote the stress 
vector over these faces by —t,;). (t(;) denotes the stress vector when +e;;) 
is the outward normal.) Then applying the principle of local equilibrium to 
the stress forces when the tetrahedron is very small we have 

Equation: 


t(,)dA = t(1)d At = t(2)dA2 = t(3)d Az 
= (tin) — taym — teaynatiayns)dA = 0 


Now let Tj; denote the qth component of tj) and ¢(,); the gth component of 
t(n) so that this equation can be written 
Equation: 


However, t/,,) is a vector and n is a unit vector quite independent of the T;; 
so that by the quotient rule the 7; are components of a second order tensor 


T. In dyadic notation we might write 
Equation: 


tin) = Ten. 


This tells us that the system of stresses in a fluid is not so complicated as to 
demand a whole table of functions t/,,) (x, n) at any given instant, but that 
it depends rather simply on n through the nine quantities T;; (x). Moreover, 
because these are components of a tensor, any equation we derive with them 
will be true under any rotation of the coordinate axis. 


Inserting the tensor expression for the stress into the momentum balance 
and using the equation of continuity and Green's theorem we have 
Equation: 
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Since all the integrals are now volume integrals, they can be combined as a 


single integrand. 
If] (eae #-¥* T)av =0 


Equation: 

However, since V is an arbitrary volume this equation is satisfied only if the 
integrand vanishes identically. 

Equation: 


> pf£+VelT 
or 
Do; 
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This is Cauchy's equation of motion and a is the acceleration. It holds for 
any continuum no matter how the stress tensor 'T is connected with the rate 
of strain. 


The symmetry of the stress tensor 


A polar fluid is one that is capable of transmitting stress couples and being 
subject to body torques, as in magnetic fluids. In case of a polar fluid we 
must introduce a body torque per unit mass in addition to the body force 
and a couple stress in addition to the normal stress t(,,.). The stress for polar 


fluids is discussed by Aris. 


A fluid is nonpolar if the torques within it arise only as the moments of 
direct forces. For the nonpolar fluid we can make the assumption either that 
angular momentum is conserved or that the stress tensor is symmetric. We 
will make the first assumption and deduce the symmetry. 


Return now to the integral linear momentum balance with the internal force 
expressed as a surface integral. If we assume that all torques arise from 
macroscopic forces, then not only linear momentum but also the angular 
momentum x x (pv) are expressible in terms of f and tn). 


Equation: 


|] fosav = |] [orav + J ws 
ele Jf focexsav-+ J [tints 


The surface integral has as its 2°” component 
Equation: 


| [extinas = Eigkh jLkpNpdS 
= [[ [esleit os) aV 


by Green's theorem. However, since 2; » = 6jp, this last integrand is 
Equation: 


Eigk(LjT pk) , = €igk&jJT php + €ijnL jr 


xx (VeT)+T, 


where T is the vector €;; T jx. 


Since v x v = 0, d(x x v)/dt = x x a, applying the transport theorem 
to the angular momentum we have 


Equation: 
=. [ [fox wav = [ [foxx aav 


Substituting back into the equation for the angular momentum and 
rearranging gives 


Equation: 
[[fr~ (pa — pf —VeT)dV = [[[t-av 
V V 


However, the left-hand side vanishes for an arbitrary volume and so 
Equation: 


T,. = 0. 


The components of T,, are (T>3 — T32), (T31 — T13), and (T12 — T>1) and 
the vanishing of these implies 
Equation: 


Ti = Ty 
so that 'T is symmetric for nonpolar fluids. 


Hydrostatic pressure 


If the stress system is such that an element of area always experiences a 
stress normal to itself and this stress is independent of orientation, the stress 
is called hydrostatic. All fluids at rest exhibit this stress behavior. It implies 
that n e T is always proportional to n and that the constant of 
proportionality is independent of n. Let us write this constant —p, then 
Equation: 


nil;; = —pnj;, hydrostatic stress 


neT = pn 


However, this equation means that any vector is a characteristic vector or 
eigenvector of T. This implies that the hydrostatic stress tensor is spherical 
or isotropic. Thus 

Equation: 


T;; = —p0;;, hydrostatic stress 
T = -pl 


for the state of hydrostatic stress. 


For a compressible fluid at rest, p may be identified with the classical 
thermodynamic pressure. On the assumption that there is local 
thermodynamic equilibrium even when the fluid is in motion this concept of 
stress may be retained. For an incompressible fluid the thermodynamic, or 
more correctly thermostatic, pressure cannot be defined except as the limit 
of pressure in a sequence of compressible fluids. We shall see later that it 
has to be taken as an independent dynamical variable. 


The stress tensor for a fluid may always be written 
Equation: 


Ty = poyt Py 
T = pl+P 


and P;; is called the viscous stress tensor. The viscous stress tensor of a 
fluid vanishes under hydrostatic conditions. 


If the external or body force is conservative (i.e., gradient of a scalar) the 
hydrostatic pressure is determined up to an arbitrary constant from the 
potential of the body force. 

Equation: 


VeT = —pf, static conditions 


f=g = gVz, uniform gravity field 
VeT = -—Vp, hydrostatic conditions 
Vp = pV®@ 
Pd 
®(p) = i = +Cy=gz+C, 


p=pgz+C;3, if p is constant 


Buoyancy (Deen, 1998) 


A consequence of the fact that pressure increases with depth in a static fluid 
is that the pressure exerts a net upward force on any submerged object. To 
calculate this force, consider an object of arbitrary shape submerged in a 
constant-density fluid as shown in the figure. The net force on this object, 
F,, is given by 

Equation: 


Pressure forces ona 
submerged object. The control 
volume in (b) is used to 
calculate the pressure force on 
the object in (a). 


where the minus sign reflects the fact that positive pressure are compressive 
(i.e., pressure acts in the —n direction). The pressure force is evaluated 
most easily by considering the situation in the figure where the solid has 
been replaced by an identical volume of fluid. Because the pressure in the 
fluid depends on depth only, this replacement of the solid does not affect the 
pressure distribution in the surrounding fluid. In particular, the pressure 
distribution on the fluid control surface in (b) must be identical to that on 
the solid surface in (a). Situation (b) has the advantage that we can apply 
the divergence theorem to the integral in the last equation, because p is a 
continuous function of position within the volume V; this is not necessarily 
true for (a), because we have said nothing about the meaning of p within a 
solid. Applying the divergence theorem and using the hydrostatic pressure 
field, we find that 

Equation: 


hy = - | fonas 
— | f [ veav = -pyev 


where pf is the (constant) density of the fluid. Thus, the upward force on 
the object due to pressure equals the weight of an equivalent volume of 
liquid; this is Archimedes’ law. An important implication of this equation 
is that F;, is independent of the absolute pressure (provided that the density 
is independent of pressure). 


The buoyancy force is the net force on the object due to the same static 
pressure variation and gravity. Evaluating the gravitational force on the 
solid body by setting p = py, the buoyancy force is found to be 
Equation: 


Fz = (9; — prev. 


Principal axes of stress and the notion of isotropy 


The diagonal terms 741, 752, 733 of the stress tensor are sometimes called 
the direct stresses and the terms 749, 751, 731, 713, T53, [32 the shear 
stresses. When there are no external or stress couples, the stress tensor is 
symmetric and we can invoke the known properties of symmetric tensors. 
In particular, there are three principal directions and referred to coordinates 
parallel to these, the shear stresses vanish. The direct stresses with these 
coordinates are called the principal stresses and the axes the principal axes 
of stress. 


An isotropic fluid is such that simple direct stress acting in it does not 
produce a shearing deformation. This is an entirely reasonable view to take 
for isotropy means that there is no internal sense of direction within the 
fluid. Another way of expressing the absence of any internally preferred 
direction is to say that the functional relation between stress and rate of 


deformation must be independent of the orientation of the coordinate 
system. We shall show in the next section that this implies that the principal 
axes of stress and rate of deformation coincide. 


The Stokesian fluid 


The constitutive equation of a non-elastic fluid satisfying the hypothesis of 
Stokes is called a Stokesian fluid. This fluid is based on the following 
assumptions. 


I. The stress tensor Tj; is a continuous function of the rate of deformation 
tensor e;; and the local thermodynamic state, but independent of other 
kinematical quantities. 

II. The fluid is homogeneous, that is, 7}; does not depend explicitly on x. 
III. The fluid is isotropic, that is, there is no preferred direction. 
IV. When there is no deformation (e;; = 0) the stress is hydrostatic, 


(Ti; = —pdiy). 


The first assumption implies that the relation between the stress and rate of 
strain is independent of the rigid body rotation of an element given by the 
antisymmetric kinematical tensor 12; ;. The thermodynamic variables, for 
example, pressure and temperature, will be carried along this discussion 
without specific mention except where it is necessary for emphasis. We are 
concerned with a homogeneous portion of fluid so the second assumption is 
that the stress tensor depends only on position through the variation of e;; 
and thermodynamic variables with position. The third assumption is that of 
isotropy and this implies that the principal directions of the two tensors 
coincide. To express this as an equation we write T;; = fi;(€pq), then if 
there is no preferred direction T;; is the same function fj; of ep, as Tj; is of 
Eng. Thus 

Equation: 


Tij = fij(€pq)- 


The fourth assumption is that the tensor 


Equation: 


Pig = Ti + pois 
vanishes when there is no motion. P;; is called the viscous stress tensor. 


Constitutive equations of the Stokesian fluid 


The arguments for the form of the constitutive equation for the Stokesian 
fluid is given by Aris and is not repeated here. The equation takes the form 
Equation: 


Ti = —pdij + Beij + Yeiren;, 
T = —pl+ Be+-yeee 


which insures that Tj; reduces to the hydrostatic form when the rate of 
deformation vanishes. 


The Newtonian fluid 


The Newtonian fluid is a linear Stokesian fluid, that is, the stress 
components depend linearly on the rates of deformation. Aris gives two 
arguments that deduce the form of the constitutive equation for a 
Newtonian fluid. 

Equation: 


Ti; = (—p al AO) 65; + 2pe;; 
where 


O=eE;=Y44,= Vev 


? 


Interpretation of the constants and yu 


Consider the shear flow given by 


Equation: 


Y= f(a2), 2 V2 = v3 = 0 


For this we have all the e;; zero except 


Equation: 

€12 = eg] = 5S" (22) = sa 
Thus 
Equation: 

Py = Px = pf’ (x2) = we 


and all other viscous stresses are zero. It is evident that pz is the 
proportionality constant relating the shear stress to the velocity gradient. 
This is the common definition of the viscosity, or more precisely the 
coefficient of shear viscosity of a fluid. 


For an incompressible, Newtonian fluid the pressure is the mean of the 
principal stresses since this is 
Equation: 


1 2 
—T;, = -p+A0+—=pO 
3 - 3h 


= —p 


For a compressible fluid we should take the pressure p as the 
thermodynamic pressure to be consistent with our ideas of equilibrium. 
Thus if we call —p the mean of the principal stresses, 

Equation: 


-(a+ ph) Vew 


2 1 Do 
2; |e 
( s") o> 


Since p, the thermodynamic pressure, is in principle known from the 
equation of state p — pis a measurable quantity. The coefficient in the 
equation is known as the coefficient of bulk viscosity. It is difficult to 
measure, however, since relatively large rates of change of density must be 
used and the assumption of linearity is then dubious. Stokes assumed that 
p = p and on this ground claimed that 


Equation: 
2 


N+ p=0 
A 


supporting this from an argument from the kinetic theory of gases. 


Assignment 5.1 
Assume a Newtonian fluid and calculate the stress tensor for the flow fields 
of assignment 4.2. Evaluate the force due to the stress on the surface y = 0. 


Equations of Motion and Energy in Cartesian Coordinates 
Topics Covered in this Chapter 


e Equations of motion of a Newtonian fluid 

e The Reynolds number 

¢ Dissipation of Energy by Viscous Forces 

¢ The energy equation 

¢ The effect of compressibility 

e Resume of the development of the equations 
e Special cases of the equations 


o Restrictions on types of motion 


= Isochoric motion 

» Irrotational motion 

= Plane flow 

» Axisymmetric flow 

= Parallel flow perpendicular to velocity gradient 


o Specialization on the equations of motion 


« Hydrostatics 

» Steady flow 

= Creeping flow 

= Inertial flow 

= Boundary layer flow 

» Lubrication and film flow 


o Specialization of the constitutive equation 


= Incompressible flow 

= Perfect (inviscid, nonconducting) fluid 
= Ideal gas 

» Pijezotropic fluid and barotropic flow 

» Newtonian fluids 


e Boundary conditions 


Surfaces of symmetry 

Periodic boundary 

Solid surfaces 

Fluid surfaces 

Boundary conditions for the potentials and vorticity 


OO > OI - O&O. 


¢ Scaling, dimensional analysis, and similarity 


o Dimensionless groups based on geometry 
o Dimensionless groups based on equations of motion and energy 


o Friction factor and drag coefficients 
e Bemoulli theorems 


o Steady, barotropic flow of an inviscid, nonconducting fluid with conservative body 
forces 

© Coriolis force 

© Jrrotational flow 

o Ideal gas 


Reading assignment 


Chapter 2&3 in BSL 
Chapter 6 in Aris 


Equations of motion of a Newtonian fluid 


We will now substitute the constitutive equation for a Newtonian fluid into Cauchy's equation 
of motion to derive the Navier-Stokes equation. 


Cauchy's equation of motion is 


Equation: 
Vi 
pou = pa = ph t+ Ty, 
or 
Dv 
= p— =pf{+VeT 
pa = pp =P ° 


The constitutive equation for a Newtonian fluid is 
Equation: 
Tiy = (p+ AO)bij + 2pess 
or 
T = (—p+A0)I+ 2ue 


The divergence of the rate of deformation tensor needs to be restated with a more meaningful 
expression. 
Equation: 


— 1 0 f dy; if Ov; 
a 2 Ox; Ox; Ox; 


1 070; re 1 0 Ov; 
2 Ox ;Ox ; 2 Ox; Ox; 


P35 1 0 
= GY Phage 
or 
he ee 
Vee = BE Nor Gg AVN) 
Thus 
Equation: 
Op fs) 
Ti,5 = Da; | (A +4) Ba, (VOY) + BV vi 
or 
Dv 
Ppp = Pf Ve + A+ H)VO + nVev 


Substituting this expression into Cauchy's equation gives the Navier-Stokes equation. 
Equation: 


Dv; Op | 0 9 
Pe TPL Be Eye os) Va 


or 


Dv 
Pa, = pf — Vp+ (A+p)VO+ pV’v 


The Navier-Stokes equation is sometimes expressed in terms of the acceleration by dividing the 
equation by the density. 


Equation: 
Dv Ov 
SS eS ES V 
oe At +(veV)v 
1 
—f a (r' + v1) VO +0V°v 


where v = p/p and A’ = A/p.v is known as the kinematic viscosity and if Stokes' relation is 
assumed X’ + v = v/3. Using the identities 
Equation: 


V-v V(Vev)-Vx(V xv) 


Vxv 


< 
| 


the last equation can be modified to give 
Equation: 


Dv 


a= — 
Dt 


=f = yee (\' + 2v) VO —vV x w. 
p 


If the body force f can be expressed as the gradient of a potential (conservative body force) and 
density is a single valued function of pressure (piezotropic), the Navier-Stokes equation can be 
expressed as follows. 


Equation: 
Dv ? 
ae ya V [Q+ P(p) — (' + 2v)O] —vV x w 
where 
Dp ! 
f=-Vi0 and P(p) = [ = 
p(p’) 


Assignment 6.1 
Do exercises 6.11.1, 6.11.2, 6.11.3, and 6.11.4 in Aris. 


The Reynolds number 


Later we will discuss the dimensionless groups resulting from the differential equations and 
boundary conditions. However, it is instructive to derive the Reynolds number Ne from the 
Navier-Stokes equation at this point. The Reynolds number is the characteristic ratio of the 
inertial and viscous forces. When it is very large the inertial terms dominate the viscous terms 
and vice versa when it is very small. Its value gives the justification for assumptions of the 
limiting cases of inviscid flow and creeping flow. 


We will consider the case of single-phase flow with conservative body forces (e.g., 
gravitational) and density a single valued function of pressure. The pressure and potential from 
the body force can be combined into a single potential. 


Equation: 
1 
f——Vp=-V2 
p 
where 
P dp 


If the change in density is small enough, the potential can be approximated by potential that has 
the units of pressure. 


Equation: 
P : : 
92 ~ —, small change in density 
p 
where 
P=p— pgz 


Suppose that the flow is characterized by a certain linear dimension, L, a velocity U, anda 
density p . For example, if we consider the steady flow past an obstacle, L may be it's diameter 
and U and p the velocity and density far from the obstacle. We can make the variables 
dimensionless with the following 

Equation: 


i FP 
Sp: 


* Vv * X * 


V — a xX ah t, P 


ace 


Vo = LV,V?=L?V? 


The conservative body force, Navier-Stokes equation is made dimensionless with these 
variables. 


Equation: 
D 
> = —-VP+(A+p)VO+pV’v 
U2 Dv" U2 * uU * uU * 
= -p—VP*+-—(A 1)VO'+-—vV"? 
PT De ep Teagan era) Ge epg ow 
L D : * kK * > * 
iaeke “+P = (A/u+1)V'e" + V2v 
bh | Dt 
Dv" * * Q_* 
Wel Ve || = Overs 
Dt 
where 
2 
Ne _ puUL pu 


wo pU/L 


The Reynolds number partitions the Navier -Stokes equation into two parts. The left side or 
inertial and potential terms, which dominates for large NRe and the right side or viscous terms, 
which dominates for small NRe. The potential gradient term could have been on the right side if 
the dimensionless pressure was defined differently, i.e., normalized with respect to (uU)/L, the 
shear stress rather than kinetic energy. Note that the left side has only first derivatives of the 
spatial variables while the right side has second derivatives. We will see later that the left side 


may dominate for flow far from solid objects but the right side becomes important in the 
vicinity of solid surfaces. 


The nature of the flow field can also be seen form the definition of the Reynolds number. The 
second expression is the ratio of the characteristic kinetic energy and the shear stress. 


The alternate form of the dimensionless Navier-Stokes equation with the other definition of 
dimensionless pressure is as follows. 
Equation: 


DV" 2k eK * > * 
Nre aa =-VP'+(A/u+1VO 4+ Vv 


P™ = P 
pu /L 


Dissipation of Energy by Viscous Forces 


If there was no dissipation of mechanical energy during fluid motion then kinetic energy and 
potential energy can be exchanged but the change in the sum of kinetic and potential energy 
would be equal to the work done to the system. However, viscous effects result in irreversible 
conversion of mechanical energy to internal energy or heat. This is known as viscous 
dissipation of energy. We will identify the components of mechanical energy in a flowing 
system before embarking on a total energy balance. 


The rate that work W is done on fluid in a material volume V with a surface S is the integral of 
the product of velocity and the force at the surface. 
Equation: 
dW 
dt 


BY b(n) dS 


§veTends 
// Ve(veT) dv 


The last integrand is rather complicated and is better treated with index notation. 
Equation: 


(viii); = Tiyvig + iT 
= Tijvij t+ vi pm = of 
= Tiyvij + see — phiri 
Ve(vet) = T: Vv + 5p —pfev 


We made use of Cauchy's equation of motion to substitute for the divergence of the stress 
tensor. The integrals can be rearranged as follows. 


Equation: 

1 1 Dv? 
al ffaere = fifa ae” 
dt 2 2 Dt 


jj paras /| Ve(veT)dv — ff fr: veav 


—p = pfev+Ve(veT)—T: Vv 


where 
T:Vv= Lis 


The left-hand term can be identified to be the rate of change of kinetic energy. The first term on 
the right-hand side is the rate of change of potential energy due to body forces. The second term 
is the rate at which surface stresses do work on the material volume. We will now focus 
attention on the last term. 


The last term is the double contracted product of the stress tensor with the velocity gradient 
tensor. Recall that the stress tensor is symmetric for a nonpolar fluid and the velocity gradient 
tensor can be split into symmetric and antisymmetric parts. The double contract product of a 
symmetric tensor with an antisymmetric tensor is zero. Thus the last term can be expressed as a 
double contracted product of the stress tensor with the rate of deformation tensor. We will use 
the expression for the stress of a Newtonian fluid. 

Equation: 


Tijvig = Ties 

—|(—p + AO) dij + 2perjlecj 
[p — A@]eis — 2wex;e:; 

= pO— 0? — 2u (oe — 26) 
= —-T:Vv 


where @ is the second invariant of the rate of deformation tensor. Thus the rate at which kinetic 
energy per unit volume changes due to the internal stresses is divided into two parts: 


i. a reversible interchange with strain energy, , 
ii. a dissipation by viscous forces, 


Equation: 
—|(A + 2p)? — 49] 


Since 9? — 2¢ is always positive, this last term is always dissipative. If Stokes' relation is used 
this term is 


Equation: 
4 
—p|—@? — 46 
: | 3 | 
for incompressible flow it is 
Equation: 
Ap®. 


(The above equation is sometimes written —4y@, where ¢ is called the dissipation function. We 
have reserved the symbol @ for the second invariant of the rate of deformation tensor, which 
however is proportional to the dissipation function for incompressible flow. Y is the symbol 
used later for the negative of the dissipation by viscous forces.) 


The energy equation 


We need the formulation of the energy equation since up to this point we have more unknowns 
than equations. In fact we have one continuity equation (involving the density and three 
velocity components), three equations of motion (involving in addition the pressure and another 
thermodynamic variable, say the temperature) giving four equations in six unknowns. We also 
have an equation of state, which in incompressible flow asserts that p is a constant reducing the 
number of unknowns to five. In the compressible case it is a relation 

Equation: 


p= f(p,T) 
which increases the number of equations to five. In either case, there remains a gap of one 
equation, which is filled by the energy equation. 


The equations of continuity and motion were derived respectively from the principles of 
conservation of mass and momentum. We now assert the first law of thermodynamics in the 


form that the increase in total energy (we shall consider only kinetic and internal energies) in a 
material volume is the sum of the heat transferred and work done on the volume. Let q denote 
the heat flux vector, then, since n is the outward normal to the surface, —q e n is the heat flux 
into the volume. Let U denote the specific internal energy, then the balance expressed by the 
first law of thermodynamics is 

Equation: 


2 
af fe (F+u)av- [[ [otevav + gveT ends ~ gqends 


This may be simplified by subtracting from it the expression we already have for the rate of 
change of kinetic energy, using Reynolds transport theorem, and Green's theorem. 


Equation: 
D 
[ [flop + Vea- ts (vv)]av =o 


Since this is valid for any arbitrary material volume, we have assumed continuity of the 
integrand 
Equation: 


We assume Fourier's law for the conduction of heat. 
Equation: 


q= —kVT 


We assume a Newtonian fluid for the dissipation of energy. 
Equation: 


T: (Vv) =—p(Vev)+Y 
Y = (A+ 2yn)0? — 4yph 


Substituting this back into the energy balance we have 
Equation: 


DU 


Physically we see that the internal energy increases with the influx of heat, the compression and 
the viscous dissipation. 


If we write the equation in the form 
Equation: 


the left-hand side can be transformed by one of the fundamental thermodynamic identities. For 
if S is the specific entropy, 
Equation: 


TdS = dU +pd(1/p) 


= dU — (p/p*)dp 


Substituting this into the last equation for internal energy gives 
Equation: 


DS 
T— = Ve(kVT)+Y. 
pT e(kVT) + 


Giving an equation for the rate of change of entropy. Dividing by 7’ and integrating over a 
material volume gives 


Equation: 
[oR — & I frsa 
= [[[|preern+ plav 
[flv (777) + =(vry + 7 |e 


= 9-S2ase |ffffcenr Fe 


The second law of thermodynamics requires that the rate of increase of entropy should be no 
less than the flux of heat divided by temperature. The above equation is consistent with this 
requirement because the volume integral on the right-hand side cannot be negative. It is zero 
only if k or VT and are zero. This equation also shows that entropy is conserved during flow if 
the thermal conductivity and viscosity are zero. 


Equation: 


DS 
— =0, when p=0, and k=0 
Dt 


Assignment 6.2 
Do exercises 6.3.1 and 6.3.2 in Aris. 


The Effect of Compressibility (Batcehlor, 1967) 


Isentropic flow. The condition of zero viscosity and thermal conductivity results in conservation 
of entropy during flow or isentropic flow. This ideal condition is useful for illustration the effect 
of compressibility on fluid dynamics. The conservation of entropy during flow implies that 
density, pressure, and temperature are changing in a reversible manner during flow. The relation 
between entropy, density, temperature, and pressure is given by thermodynamics. 


Equation: 
Os Os 
= | —] dT —]d 
= ey (Fr) @ 


C. 
SSE ares (52) dp 
Pp 


T OT 
Cp B 
= —dT——d 
T ‘p 
where 


These relations may be combined with the condition that the material derivative of entropy is 
zero to obtain a relation between temperature and pressure during flow. 
Equation: 


o DT _ BT Dp 
ED DE 


, when w = 0, and k=0 


The equation of state expresses the density as a function of temperature and pressure. During 
isentropic flow the pressure and temperature are not independent but are constrained by 
constant entropy or adiabatic compression and expansion. The density in this case is given as 
Equation: 


p = p(p, 8) 


We now have as many equations as unknowns and the system can be determined. The 
simplifying feature of isentropic flow is that exchanges between the internal energy and other 
forms of energy are reversible, and internal energy and temperature play passive roles, merely 
changing in response to the compression of a material element. The continuity equation and 
equation of motion governing isotropic flow may now be expressed as follows. 

Equation: 


1 Dp 

pe ee 
Dv 
— =pf-—V 
Pay p 

where 


The physical significance of the parameter c, which has the dimensions of velocity, may be seen 
in the following way. Suppose that a mass of fluid of uniform density pz, is initially at rest, in 
equilibrium, so that the pressure p, is given by 

Equation: 


Vio = Pof. 


The fluid is then disturbed slightly (all changes being isentropic), by some material being 
compressed and their density changed by small amounts, and is subsequently allowed to return 
freely to equilibrium and to oscillate about it. (The fluid is elastic, and so no energy is 
dissipated, so oscillations about the equilibrium are to be expected.) The perturbation quantities 
P1 (= Pp — Po) and pi (= p — po) and v are all small in magnitude and a consistent 
approximation to the continuity equation and equations of motion is 


Equation: 
1 Op; 
—— —+p.Vev=0 
c2 Ot e 
Ov 
o— = pif — V 
p By Pl Pi 


where Cz is the value of c at p = p. . On eliminating v we have 
Equation: 


1 0° p1 


3 Bt? = V* pi = Vie (pif) 


The body forces commonly arise from the earth's gravitational field, in which case the 
divergence is zero and the last term is negligible except in the unlikely event of a length scale of 
the pressure variation not being small compared with c2/g (which is about 1.2 x 104m for air 
under normal conditions and is even larger for water). Thus under these conditions the above 
equation reduces to the wave equation for p; and p, satisfies the same equation. The solutions 
of this equation represents plane compression waves, which propagate with velocity c, and in 
which the fluid velocity v is parallel to the direction of propagation. In another words, c, is the 
speed of propagation of sound waves in a fluid whose undisturbed density is po. 


Conditions for the velocity distribution to be approximately solenoidal. The assumption of 
solenoidal or incompressible fluid flow is often made without a rigorous justification for the 
assumption. We will now reexamine this assumption and make use of the results of the previous 
section to express the conditions for solenoidal flow in terms of identifiable dimensionless 


groups. 


The condition of solenoidal flow corresponds to the divergence of the velocity field vanishing 
everywhere. We need to characterize the flow field by a characteristic value of the change in 
velocity U with respect to both position and time and a characteristic length scale over which 
the velocity changes L. The spatial derivatives of the velocity then is of the order of U/L. The 
velocity distribution can be said to be approximately solenoidal if 


Equation: 
U 
|\Vevl< E 
ie., if 
1 Dp U 
pDt ~ L£ 


For a homogeneous fluid we may choose p and the entropy per unit mass S as the two 
independent parameters of state, in which case the rate of change of pressure experienced by a 
material element can be expressed as 


Equation: 
p = v(p, 8) 
Dp _ ,Dp(Op\ DS 
Dit ° Dt\ dS), Dt’ 


The condition that the velocity field should be approximately solenoidal is 
Equation: 


1 Dp 1 (3) DS U 


< : 
pe? Dt pe? \ OS} , Dt L 


This condition will normally be satisfied only if each of the two terms on the left-hand side has 
a magnitude small compared with U/L. We will now examine each of these terms. 


I. When the condition 
Equation: 


is satisfied, the changes in density of a material element due to pressure variations are 
negligible, i.e., the fluid is behaving as if it were incompressible. This is by far the more 
practically important of the two requirements for v to be a solenoidal vector field. In 
estimating | Dp/Dt| we shall lose little generality by assuming the flow to be isentropic, 
because the effects of viscosity and thermal conductivity are normally to modify the 
distribution of pressure rather than to control the magnitude of pressure variation. We may 
then rewrite the last equation with the aid of equations of motion of an isentropic fluid 
derived in the last section. 


Equation: 
Dp Op 
= V; 
Dt Ot ee ee 
op ODE fas cll ee 
= Ot at 
= OP 2: Pepe 
~ aE? 20at 
Thus 
Equation: 


1 Op , vef 1 Dv? U 
pee Ot. 2c? Dt iy 


Showing that in general three separate conditions, viz. that each term on the left-hand side 
should have a magnitude small compared with U/L if the flow field is to be 
incompressible. 


° I (i) Consider first the last term on the left-hand side of the above equation. The order 
of magnitude of Dv? / Dt will be the same as that of 0v?0/t or v e Vv" (i.e., U?/L), 
which ever is greater. Thus the condition arising from this term can be expressed as 
Equation: 


IT: 


o I (ii) The magnitude of the partial derivative of pressure with respect to time depends 
directly on the unsteadiness of the flow. Let us suppose that the flow field is 
oscillatory and that v is a measure of the dominant frequency. The rate of change of 
momentum is then the order of pUv . Since the pressure gradient is the order of the 
rate of change of momentum, the spatial pressure variation over a region of length L 
is pLUV . Since the pressure is also oscillating, the magnitude of Op/Ot is then 
OLUv”. Thus the condition that the first term be small compared to U/L is 
Equation: 


all Bg 
Pe 


ily 


This condition is equivalent to the condition that the length of the system should be 
small enough that a pressure transient due to compression is felt instantaneously 
throughout the system. 

o (iii) If we regard the body forces arising from gravity, the term from the body forces, 
v ef/c?, has a magnitude of order gU/c?, so the condition that it be small compared 
toU/Lis 
Equation: 


el? 


<i 
C2 


This condition is equivalent to the condition that the length of the system should be 
small enough that a pressure transient due to compression is felt instantaneously 
throughout the system. 


This shows that the condition is satisfied provided the difference between the static-fluid 
pressure at two points at vertical distance L apart is a small fraction of the absolute 
pressure, i.e., provided the length scale L characteristic of the velocity distribution is small 
compared to p/pg, the 'scale height’ of the atmosphere, which is about 8.4 km for air under 
normal conditions. The fluid will thus behave as if it were incompressible when the three 
conditions I(i), I(ii), and I(iii) are satisfied. The first is not satisfied in near sonic or 
hypersonic gas dynamics, the second is not satisfied in acoustics, and the third is not 
satisfied dynamical meteorology. 

The second condition necessary for incompressible flow is that arising from entropy. This 
condition requires that variation of density of a material element due to internal dissipative 
heating or due to molecular conduction of heat into the element be small. We will show 


later how the small variation of density leading to natural convection can be allowed by yet 
assume incompressible flow. 


Resume of the development of the equations 


We have now obtained a sufficient number of equations to match the number of unknown 
quantities in the flow of a fluid. This does not mean that we can solve them nor even that the 
solution will exists, but it certainly a necessary beginning. It will be well to review the 
principles that have been used and the assumptions that have been made. 


The foundation of the study of fluid motion lies in kinematics, the analysis of motion and 
deformations without reference to the forces that are brought into play. To this we added the 
concept of mass and the principle of the conservation of mass, which leads to the equation of 
continuity, 

Equation: 


<P 5 p(Vev) = 2 + Ve (pv) =0 


An analysis of the nature of stress allows us to set up a stress tensor, which together with the 
principle of conservation of linear momentum gives the equations of motion 
Equation: 


Dv 
—=-=pf+VeT. 
Ppp re 


If the conservation of moment of momentum is assumed, it follows that the stress tensor is 
symmetric, but it is equally permissible to hypothesize the symmetry of the stress tensor and 
deduce the conservation moment of momentum. For a certain class of fluids however (hereafter 
called polar fluids) the stress tensor is not symmetric and there may be an internal angular 
momentum as well as the external moment of momentum. 


As yet nothing has been said as to the constitution of the fluid and certain assumptions have to 
be made as to its behavior. In particular we have noticed that the hypothesis of Stokes that leads 
to the constitutive equation of a Stokesian fluid (not-elastic) and the linear Stokesian fluid 
which is the Newtonian fluid. 

Equation: 


Ti; = (—p + a) dj; + Bei + yeerj,  Stokesian fluid 
Ti; = (—- 


(—p + AO) dij + 2pyei;, Newtonian fluid. 


The coefficients in these equations are functions only of the invariants of the rate of 
deformation tensor and of the thermodynamic state. The latter may be specified by two 
thermodynamic variables and the nature of the fluid is involved in the equation of state, of 
which one form is 


Equation: 


p= f(p,T). 


If we substitute the constitutive equation of a Newtonian fluid into the equations of motion, we 
have the Navier-Stokes equation. 
Equation: 


D 
Pa, = pf —Vp+ (A+ n)V(Vev) + pov 


Finally, the principle of the conservation of energy is used to give an energy equation. In this, 
certain assumptions have to be made as to the energy transfer and we have only considered the 
conduction of heat, giving 

Equation: 


D 
= = Ve(kVT)—p(Vev)+Y 


These equations are both too general and too special. They are too general in the sense that they 
have to be simplified still further before any large body of results can emerge. They are too 
special in the sense that we have made some rather restrictive assumptions on the way, 
excluding for example elastic and electromagnetic effects. 


Special cases of the equations 


The full equations may be specialized is several ways, of which we shall consider the 
following: 


i. restrictions on the type of motions, 
ii. specializations on the equations of motion, 
iii. specializations of the constitutive equation or equation of state. 


This classification is not the only one and the classes will be seen to overlap. We shall give a 
selection of examples and of the resulting equations, but the list is by no means exhaustive. 


Under the first heading we have any of the specializations of the velocity as a vector field. 
These are essentially kinematic restrictions. 


e (ia) Isochoric motion. (i.e., constant density) The velocity field is solenoidal 
Equation: 


The equation of continuity now gives 
Equation: 
Dp | 
Dt 


0 


that is, the density does not change following the motion. This does not mean that it is 
uniform, though, if the fluid is incompressible, the motion is isochoric. The other equations 
simplify in this case for we have a , 8 , and y of the constitutive equations functions of 
only and W of the invariants of the rate of deformation tensor. In particular for a 
Newtonian fluid 

Equation: 


Li = po Shey 
Dv 


= Spf ="V; Vv? 
p Di p pt+ pV'v 
The energy equation is 
Equation: 
DU 
—— =Ve(kVT)+Y, 
Pe (kVT) 


and for a Newtonian fluid 
Equation: 


Y = —4po. 


Because the velocity field is solenoidal, the velocity can be expressed as the curl of a 
vector potential. 
Equation: 


v=VxA. 


The Laplacian of the vector potential can be expressed in terms of the vorticity. 
Equation: 
w= Vxv 
= Vx(Vx A) 
= V(VeA)—-V’A 
—V°A, if VeA=0 
If the body force is conservative, i.e., gradient of a scalar, the body force and pressure can 


be eliminated from the Navier-Stokes equation by taking the curl of the equation. 
Equation: 


D 
eM Ae Vee 
Dt 


where v is the kinematic viscosity. 

Isochoric motion is a restriction that has to be justified. Because it is justified in so many 
cases, it is easier to identify the cases when it does not apply. We showed during the 
discussion of the effects of compressibility that compressibility or non-isochoric is 
important in the cases of significant Mach number, high frequency oscillations such as in 
acoustics, large dimensions such as in meteorology, and motions with significant viscous 
or compressive heating. 

(ib) Irrotational motion. The velocity field is irrotational 

Equation: 


w=Vxv=0 


It follows that there exists a velocity potential (x,t) from which the velocity can be derived 
as 
Equation: 


v=V¢ 


and in place of the three components of velocity we seek only one scalar function. (Note 
that some authors express the velocity as the gradient of a scalar and others as the negative 
of a gradient of a scalar. We will use either to conform to the book from which it was 
extracted.) The continuity equation becomes 


Equation: 

Dp 

—+pV'd=0 

pe teY 
so that for an isochoric (or incompressible), irrotational motion, is a potential function 
satisfying 
Equation: 

V’°d=0 
The Navier-Stokes equations become 
Equation: 
pee + (79) iS eee (N' + 2v)V(V7¢) 
Oot 2 p 

In the case of an irrotational body force f = —V2 and when p is a function only of , this 


has an immediate first integral since every term is a gradient. Thus if P (p) = f dp/p, 
Equation: 


a6 1 


5 3 (Vb) +2+ P(p) — (0+ 2) Vb = (t) 


is a function of time only. 

Irrotational motions with finite viscosity are only very special motions because the no-slip 
boundary conditions on solid surfaces usually will cause generation of rotation. Usually 
irrotational motion is associated with inviscid fluids because the no-slip boundary 
condition then will not apply and initially irrotational motion will remain irrotational. 

(ic) Complex lamellar motions, Betrami motions, ect. These names can be applied when 
the velocity field is of this type. Various simplifications are possible by expressing the 
velocity in terms of scalar fields. We shall not discuss them further here. 

(id) Plane flow. Here the motion is restricted to two dimensions which may be taken to be 
the 012 plane. Then v3 = 0 and x3 does not occur in the equations. Also, the vector 
potential and the vorticity have only one nonzero component. 

Incompressible plane flow. Since the flow is solenoidal, the velocity can be expressed as 
the curl of the vector potential. The nonzero component of the vector potential is the 
stream function. 

Equation: 


Vv VxA 
Ue => Cpa =e 
= Ow Op 
(v1, V2, v3) = (=. Be? 0) 


The vorticity has only a single component, that in the 03 direction, which we will write 
without suffix 


Equation: 
w = Vxv 
Wi = EijkURj, J,k=1,2 
~*2 Ov Ov1 
w = Ox, _ Ox2 
= —Vp 


If the body force is conservative, i.e., gradient of a scalar, then the body force and pressure 
disappear from the Navier-Stokes equation upon taking the curl of the equations. In plane 
flow 

Equation: 


Thus for incompressible, plane flow with conservative body forces, the continuity equation 
and equations of motion reduce to two scalar equations. 

Incompressible, irrotational plane motion. A vector field that is irrotational can be 
expressed as the gradient of a scalar. Since the flow is incompressible, the velocity vector 
field is solenoidal and the Laplacian of the scalar is zero, i.e., it is harmonic or an 
analytical function. 

Equation: 


v= Vo 
Vev 0 
Vo 


Since the flow is incompressible, it also can be expressed as the curl of the vector 
potential, or in plane flow as derivatives of the stream function as above. Since the flow is 
irrotational, the vorticity is zero and the stream function is also an analytical function, i.e., 
v=VxA, 0=w=Vxv=-VA4V(VeA),> V4 =0. Thus 
Equation: 


— Ob _ oy 
ss 1 = Ox 1 = Ox 2 

meh. SEND 9 ote Oy 
u 2 = Ox 2 = Ox 1 


These relations are the Cauchy-Riemann relations show that the complex function 

f = ¢@+71yis an analytical function of z = x; + 722. The whole resources of the theory 
of functions of a complex variable are thus available and many solutions are known. 
Steady, plane flow. If the fluid is compressible but the flow is steady (i.e., no quantity 
depends on t) the equation of continuity becomes 

Equation: 


Movi)», Cpe) 
Ox, 0x2 


A stream function can again be introduced, this time in the form 
Equation: 


2, AO _1 Oy 
p Ox,’ ” p Ox 


Cal 


The vorticity is now given by 
Equation: 


(ie) Axisymmetric flows. Here the flow has an axis of symmetry such that the flow field 
can be expressed as a function of only two coordinates by using curvilinear coordinates. 
The curl of the vector potential and velocity has only one non-zero component and a 
stream function can be found. 

(if) Parallel-flow perpendicular to velocity gradient. If the flow is parallel, i.e., the 
streamlines are parallel and are perpendicular to the velocity gradient, then the equations 
of motion become linear in velocity if the fluid is Newtonian. 

Equation: 


If veVv=0,_ then by & and 


Dt 
p& =f+VeT 
The second type of specializations are limiting cases of the equations of motion. 


(iia) Hydrostatics. When there is no flow, the only non-zero terms in the equations of 
motion are the body forces and pressure gradient. 
Equation: 


pi=Vp 


If the body force is conservative, i.e., the gradient of a scalar, then the hydrostatic pressure 
can be determined from this scalar. 
Equation: 

F=V 0D 

Vp = pV2Q 

V(é— 2) =0 

@ = (92+ constant 

where 


Pp 
$=f2 


If the body force is due to gravity then 
Equation: 


Q= 9z 


where 2 is the elevation above a datum such as the mean sea level. 

(iib) Steady flow. Examples of this have already been given and indeed it might have been 
considered as a restriction of the first class. All partial derivatives with respect to time 
vanish and the material derivative reduce to the following 

Equation: 


D 
rs aaa 
In particular, the continuity equation is 
Equation: 
Ve(pv)=0 


so that the mass flux field is solenoidal. 
(iic) Creeping flow. It is sometimes justifiable to assume that the velocity is so small that 
the square of velocity is negligible by comparison with the velocity itself. This linearizes 


the equations and allows them to be solved more readily. For example, the Navier-Stokes 
equation becomes 
Equation: 


Ov 


Pans f—Vp+ (A+ p)V(Vev) +pV?v 


In particular, for steady, incompressible, creeping flow with conservative body forces 
Equation: 


VPS Vey. 
where 
P=p-—pgz 


However, since the continuity equation is V e v = 0, we have 
Equation: 


V-P=0 


or P is a harmonic potential function. This is the starting point for Stokes’ solution of the 
creeping flow about a sphere and for its various improvements. 

One may ask, "How small must the velocity be in order to neglect the nonlinear terms?" To 
answer this question, we need to examine the value of the Reynolds number. However, this 
time the pressure and body forces will be made dimensionless with respect to the shear 
forces rather than the kinetic energy. 

Equation: 


* 
Veo — 
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ae Pa, VPS 
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V=LV, Vo SPW 


The dimensionless Navier-Stokes equation for incompressible flow is now as follows 
Equation: 


* KK x2 * 
Nr? = -VP™4+Vv 
where 
= pUL _ pU? 
Nre — yw  pu/o 


Creeping flow is justified if the Reynolds number is small enough to neglect the left-hand 
side of the above equation. If the dimensionless variables and their derivatives are the 
order of unity, then creeping flow is justified if the Reynolds number is small compared to 
unity. 

Another specialization of the equations of motion where the equations are made linear 
arises in stability theory when the basic flow is known but perturbed by a small amount. 
Here it is the squares and products of small perturbations that are regarded as negligible. 


¢ (iid) Inertial flow. The flow is said to be inviscid when the inertial terms are dominant and 
the terms with viscosity in the equations of motion can be neglected. We can examine the 
conditions when this may be justified from the dimensionless equations of motion with the 
pressure and body forces normalized with respect to the kinetic energy. 
Equation: 


NRe bra + vP*| = (A/u+1) V0" + Vv" 
where 
_ pUL _ pU? 
INES? = wo pU/L 
P= mie 


The limit of inviscid flow may occur when the Reynolds number is becomes very large 
such that the right-hand side of the above equation is negligible. Notice that if the right- 
hand side of the equation vanishes, the equation goes from being second order in spatial 
derivatives to first order. The differential equation goes from being parabolic to being 
hyperbolic and the number of possible boundary condition decreases. The no-slip 
boundary condition at solid surfaces can no longer apply for inviscid flow. These types of 
problems are known as singular perturbation problems where the differential equations are 
first order except near boundaries where they become second order. Physically, the form 
drag dominates the skin friction in inertial flow. The macroscopic momentum balance is 
described by Bernoulli theorems. 
Notice that inviscid flow is necessary for irrotational flow past solid objects but inviscid 
flow may be rotational. If fact much of the classical fluid dynamics of vorticity is based on 
inviscid flow. 

¢ (iie) Boundary-layer flows. Ideal fluid or inviscid flow may be assumed far from an object 
but real fluids have no-slip boundary conditions on solid surfaces. The result is a 
boundary-layer of viscous flow with a large vorticity in a thin layer near solid surfaces that 
merges into the ideal fluid flow at some distance from the solid surface. It is possible to 
neglect certain terms of the equations of motion compared to others. The basic case of 
steady incompressible flow in two dimensions will be outlined. If a rigid barrier extends 
along the positive 01 axis the velocity components v; and v2 are both zero there. In the 
region distant from the axis the flow is v1 = U (x1), v2 = 0, and may be expected to be of 
the form shown in Fig. 6.1, in which v1 differs from U and v2 from zero only within a 
comparatively short distance from the plate. To express this we suppose L is a typical 
dimension along the plate and 6 a typical dimension of this boundary layer and that V; and 
V2 are typical velocities of the order of magnitude of v1 and v2. We then introduce 
dimensionless variables 
Equation: 

* L1 * Bip) * U1 * v2 


4 = ad aa ae — eee, | 


which will be of the order of unity. This in effect is a stretching upward of the coordinates 
so that we can compare orders of magnitude of the various terms in the equations of 
motion, for now all dimensionless quantities will be of order of unity. It is assumed that 


0 << L, but the circumstances under which this is valid will become apparent later. It is 
also assumed that the functions are reasonably smooth and no vast variations of gradient 
occur. The equation of continuity becomes 

Equation: 


Vi du" Vo dv" 


L éx* 6 dy 


which would lose its meaning if one of these terms were completely negligible in 
comparison with the other. It follows that 


Equation: 
é 
Vo = (vi z) 


where the symbol O means "is of the order of." The Navier-Stokes equations become 
Equation: 


Vie * Ou Vivo .* Ou _ Po Op at veal O7u a L? Gu 
L Ox* 6 oy" pL dx" LT? \ dx"? 5° Oy” 
and 
Vivo, * dv , Ve * du _ po Op , 1 Va ( & Oy" | O?u" 
L en) dy pd Oy 1 BNL? Oe? ' Oy? 


In the first of these equations the last term on the right-hand side dominates the Laplacian 
and 0?u' /Ox’? can be neglected. Dividing through by V;2/L we see that the other terms 
will be of the same order of magnitude provided 

Equation: 


Lv 


= = OL 
py, OW) 


Po=pV, and 


Inserting these orders of magnitude in the second equation and the condition that 6 << D 
results in the pressure gradient Op/Oy being the dominant term. Thus p is a function of «x 
only. Returning to the original variables, we have the equations 


Equation: 
Ov} Ov2 = 
Ox, a Ox ~— 0 
Ovi | Ov, __ 1 dp | 07v1 
UL dx, | v2 Or. p dx | y da? 


The circumstances under which these simplified equations are valid are given by the term 
above that we assumed to be of order of unity, which we rewrite as 
Equation: 


Since it was assumed that 6 << L this equation shows that this will be the case if 

v < V\L. The assumption that the pressure is O(pVj2) is consistent with the assumption 
that the outer flow is inertial denominated flow. 

(iif) Lubrication and film flow. Lubrication and film flow is another case where one 
dimension is small compared to the other dimensions. Lubrication flow is usually between 
two solid surfaces with the relative velocity between the two surfaces specified. Film flow 
is usually between a solid and fluid surfaces or between two fluid interfaces and the 
driving force for flow is usually buoyancy or the Laplace pressure between curved 
interfaces. Because of the one dimension being small, lubrication and film flow is about 
always at low Reynolds numbers. 

The equations of motion are specialized by normalizing the variables with characteristic 
quantities of the flow and dropping the terms that are negligible. Incompressible flow with 
low Reynolds number of a Newtonian fluid is assumed. Let the direction normal to the 
film be the 03 direction. The double subscript, 12, will be used to denote the two 
coordinate directions in the plane of the film. Let h, be a characteristic film thickness and 
L be a characteristic dimension in the plane of the film and h,/L << 1. The variables are 
normalized to be order of unity. 


Equation: 
* * 
— £12 — «x3 eee Pint 
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where 


P,o=ApglL or Po=% or P, = £(42 L) 


The film boundaries are material surfaces and the velocity perpendicular to the plane of the 
film is the rate of change of film thickness. 
Equation: 

0 
v3|, = oF 
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Substitution into the integral of the equation of continuity result in the following. 
Equation: 


The characteristic time can be chosen as to make the dimensionless group equal to unity. 


Equation: 


The dimensionless variables can now be substituted into the equations of motion for zero 
Reynolds number with gravitational body force and Newtonian fluid. For the components 
in the plane of the film, 

Equation: 


+ ot (#U] ont , [aut Ov 
P. (2) h2 Ox a 
The dimensionless group in the last term can be made equal to unity because of the two 


characteristic quantities, P, and U, only one is specified and the other can be determined 
from the first. 


Equation: 
BUL _ 
Phe 1 
and 
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The equations of motion in the plane of the film is now 
Equation: 


The dimensionless variables can now be substituted into the equation of motion 
perpendicular to the plane of the film. 
Equation: 


0 = oP" ipa *V2u" 4 UL | [ow 12 8%, 
iit Pr | | L 12U3 L 
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This last equation shows that the pressure is approximately uniform over the thickness of 
* 


the film. Thus the velocity profile of vj. can be determined by integrating the equation of 
motion in the plane of the film twice and applying the appropriate boundary conditions. 
The average velocity in the film can be determined by integrating the velocity profile 
across the film. 

Certain of the specializations based on the constitutive equation or equation of state have 
turned up already in the previous cases. We mention here a few important cases. 


¢ (iiia) Incompressible fluid. An incompressible fluid is always isochoric and the 
considerations of (ia) apply. It should be remembered that for an incompressible fluid the 
pressure is not defined thermodynamically, but is an variable of the motion. 

e (iiib) Perfect fluid. A perfect fluid has no viscosity so that 
Equation: 


r= =pl 


ptr =pf—Vp 
If, in addition, the fluid has zero conductivity the energy equation becomes 
Equation: 
DS 
pas a) 
Dt 
and the flow is isentropic. 


¢ (iiic) Ideal gas. An ideal gas is a fluid with the equation of state 
Equation: 


p=pRT 


The entropy of an ideal gas is given by 


Equation: 
dT 
= 1 Cy — —-RIin p 
t 


which for constant specific heats gives 
Equation: 


Pp — eS/ce p” 


¢ (iiid) Piezotropic fluid and barotropic flow. When the pressure and density are directly 
related, the fluid is said to be piezotropic. A simple relation between p and allows us to 
write 
Equation: 


Pp 
<Vp=VP(o)=Vv f @ 
p p 


e (iiie) Newtonian fluids. Here the assumption of a linear relation between stress and strain 
leads to the constitutive equation 
Equation: 


T = (—p+A0O)I+2ue 


The equations of motion becomes the Navier-Stokes equations. 


Assignment 6.3 
For 2-D, incompressible flow, prove that 7 = [nevds 


Assignment 6.4 
Carry out the steps in specializing the continuity and equations of motion for boundary layer 
and lubrication or film flows. 


Assignment 6.5 
Derive the equations for the propagation of acoustic waves. 


Assignment 6.6 
Derive the simplifications that are possible to the continuity equation and equations of motion 
of Newtonian fluids for the following cases: 


1. Incompressible flow 

2. Irrotational flow 

3. Incompressible and irrotational flow 

4. Low Reynolds number. 

5. One dimension small compared to other dimensions. 
6. Laminar boundary layer. 


Boundary conditions 


The flow field is often desired for a finite region of space that is bounded by a surface. 
Boundary conditions are needed on these surfaces and at internal interfaces for the flow field to 
be determined. The boundary conditions for temperature and heat flux are continuity of both 
across internal interfaces that are not sources or sinks and either a specified temperature, heat 
flux, or a combination of both at external boundaries. 


Surfaces of symmetry. Surfaces of symmetry corresponds to reflection boundary conditions 
where the normal component of the gradient of the dependent variables are zero. Thus surfaces 
of symmetry have zero momentum flux, zero heat flux, and zero mass flux. Because the 
momentum flux is zero, the shear stress is zero across a surface of symmetry. 


Periodic boundary. Periodic boundaries are boundaries where the dependent variables and its 
derivatives repeat themselves on opposite boundaries. The boundaries may or may not be 
symmetry boundary conditions. An example of when periodic boundaries are not symmetry 
boundaries are the boundaries of 9 = 0 and 8 = 27 of a non-symmetric system with cylindrical 
polar coordinate system. 


Solid surfaces. A solid surface is a material surface and kinematics require that the mass flux 
across the surface to be zero. This requires the normal component of the fluid velocity to be that 
of the solid. The tangential component of velocity depends on the assumption made about the 
fluid viscosity. If the fluid is assumed to have zero viscosity the order of the equations of 
motion reduce to first order and the tangential components of velocity can not be specified. 


Viscous fluids stick to solid surfaces and the tangential components of velocity is equal to that 
of the solid. Exception to the 'no-slip' boundary conditions is when the mean free path of as gas 
is similar to the dimensions of the solid. An example is the flow of gas through a fine pore 
porous media. 


Porous surface. A porous surface may not be a no-flow boundary. Flux through a porous 
material is generally described by Darcy's law. 


Fluid surfaces. If there is no mass transfer across a fluid-fluid interface, the interface is a 
material surface and the normal component of velocity on either side of the interface is equal to 
the normal component of the velocity of the interface. The tangential component of velocity at a 
fluid interface is not known apriori unless the interface is assumed to be immobile as a result of 
adsorbed materials. The boundary condition at fluid interfaces is usually jump conditions on the 
normal and tangential components of the stress tensor. Aris give a thorough discussion on the 
dynamical connection between the surface and its surroundings. If we assume that the 
interfacial tension is constant and that it is possible to neglect the surface density and the 
coefficients of dilational and shear surface viscosity then the jump condition across a fluid-fluid 
interface is 

Equation: 


[Tis] ni = —2H on; 
[Tlen=-2Hon 


where the bracket denotes the jump condition across the interface, H is the mean curvature of 
the interface, and a is the interfacial or surface tension. The jump condition on the normal 
component of the stress is the jump in pressure across a curved interface given by the Laplace- 
Young equation. The tangential components of stress are continuous if there are no surface 
tension gradients and surface viscosity. Thus the tangential stress at the clean interface with an 
inviscid fluid is zero. 


For boundary conditions at a fluid interface with adsorbed materials and thus having interfacial 
tension gradients and surface viscosity, see Chapter 10 of Aris and the thesis of Singh (1996). 


Boundary_conditions for the potentials and vorticity. Some fluid flow problems are more 
conveniently calculated through the scalar and vector potentials and the vorticity. 
Equation: 


The boundary condition on the scalar potential is that the normal derivative is equal to the 
normal component of velocity. 
Equation: 


neV¢d = oe 


The boundary condition on the vector potential is that the tangential components vanish and the 
normal derivative of the normal component vanish (Hirasaki and Hellums, 1970). Wong and 
Reizes (1984) introduced a method where the need for calculation of the scalar potential is 
replaced by the use of an irrotational component of velocity. 

Equation: 


A(t) = 0 


In two dimensional or axisymmetric incompressible flow, it is not necessary to have a scalar 
potential and the single nonzero component of the vector potential is the stream function. The 
boundary condition on the stream function for flow in the x1, x2 plane of Cartesian coordinates 
is 


Equation: 
nev = neVxA 
_ Ow ow 
= n1 Bay ng Oa. 
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The boundary condition on the normal component of the vorticity of a fluid with a finite 
viscosity on a solid surface is determined from the tangential components of the velocity of the 
solid. It is zero if the solid is not rotating. If the boundary is an interface between two viscous 
fluids then the normal component of vorticity is continuous across the interface (C. Truesdell, 
1960). If the interface is with an inviscid fluid, the tangential components of vorticity vanishes 
and the normal derivative of the normal component vanishes for a plane interface (Hirasaki, 
1967). 


If the boundary is at along a region of space for which the velocity field is known, the vorticity 
can be calculated from the derivatives of the velocity field. 


If the boundary is a surface of symmetry, the tangential components of vorticity must vanish 
because the normal component of velocity and the normal derivative of velocity vanish. The 
normal derivative of the normal component of vorticity vanishes from the solenoidal property 
of vorticity. 


Scaling, Dimensional Analysis, and Similarity 


We have already seen some examples of scaling and dimensional analysis when we determined 
when the continuity equations and equations of motion could be simplified. The concept of 
similarity states that the solution of transport problems do not need to be determined separately 
for each value of the parameters. Rather the variables and parameters can be grouped into 
dimensionless variables and dimensionless numbers and the solution will have fewer degrees of 
freedom. Also, in some cases the partial differential equations can have the independent 
variables combined to fewer independent variables and be expressed as ordinary differential 
equations. The concept of similarity does not apply only to mathematical solutions but is also 
used to design physical analogs of systems on a smaller scale or with different transport 
mechanism. For example, before numerical simulation the streamlines and pressure gradients 
for flow in petroleum reservoirs were studied by electrical conduction on a laboratory scale 
model that is geometrically similar. 


Dimensionless groups based on geometry. The aspect ratio is the ratio of the characteristic 
lengths of the system. The symbol a is normally used to denote the aspect ratio. We saw how 
the aspect ratio simplified the equations of motion for boundary layer flow and lubrication or 
film flow. The table below lists some examples of aspect ratio expressions for different 
problems. 


Examples of aspect ratio, a 


6/L boundary layer flow 

ho/L lubrication or film flow 

pion ee width/length 

D/L entrance effects in pipe flow 
Dmaz/Dmin eccentricity ratio 


It may be difficult to define the characteristic length of an irregularly shaped conduit or object. 
The characteristic dimension for an irregular conduit or object can be determined by the 
hydraulic diameter. 

Equation: 
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+, 2 —D object 


vee 3 — D object 


where A, and P,, are the cross-sectional area and wetted perimeter of the conduit and V and 
Aw are the volume and wetted area of an object. The definitions reduce to the diameter of a 
cylinder or sphere for regular objects. (The reader should be aware that some definitions such as 
the hydraulic radius in BSL may not reduce to the dimension of the regular object.) The 
hydraulic diameter may provide length scales but exact similarity is not satisfied unless the 
conduits or objects are geometrically similar. 


Dimensionless groups based on equations of motion and energy. We derived earlier the 
Reynolds number from the equations of motion and dimensionless groups from the energy 
equation when compressibility is important. We will discuss these further and additional 
dimensionless groups. 


Interpretation of the Reynolds number 


puUL L . . *,e 
. in basic definition 
pu? kinetic energy 
wU/L shear stress 
~ lpve(Vv)| inertial force 
|pV2v| viscous force 


The Reynolds number can be interpreted as a ratio of kinetic energy to shear stress. We will see 
later that some of the dimensionless numbers differ by whether it is normalized with respect to 
the kinetic energy or the shear stress term, i.e., it is a product of the Reynolds number and 
another dimensionless number. 


Suppose the characteristic value of the body force is pgL and the characteristic value of 


pressure is o/L, i.e., due to capillary forces. The dimensionless Navier-Stokes equation can be 
expressed as follows. 


Equation: 
Dv* gL |/(f o 1 2 
= V*p*+|——| Vv 
pe = [oe |(5) ~ [pam]? Le] 7 


We can now define dimensionless groups that include gravity or buoyancy forces and gravity 
forces. 


Dimensionless groups based on gravity and capillarity 


2 U2 
Nr, = a = & Froude number 


ae aE Weber Number 


U U/L 
Nog=# _— wU/L _ Nwe 


o Gis Nas capillary number 
pgl? _—pglh_ __ Nre ; 
pU  pwU/L Np, gravity number 


_ pgl? _ pol Nwe 
Nao = = = SE ln ae Bond number 


We derived scale factors that have to be small in order to neglect the effects of compressibility. 
They are summarized here. 


Dimensionless groups necessary for incompressible flow 


Nua = 4% Mach number 

a (frequency length)/sonic velocity 
oe eee density change due to body f 
Ga ag ensity change due to body forces 


Friction factor and drag coefficients. Friction factor and drag coefficients are the force on the 
wall of a conduit or on an object normalized with respect to kinetic energy. There are some 
ambiguities in the literature that one should be aware of. 


Friction factors 


Tw 
fsp = pU2 Stanton-Pannel 


Fanning 


— _2Tw 
fr= pU;, 


_— _8tw ‘ 
fow = paz Darcy-Weisbach 
2D AP/L 
fMoody=Gaz = 4fr = fow Moody 
_ F frictiont+ Farag ae 
farag = A, (1/2pUz) drag coefficient 


F rictiont Fg a. 
C, eed rag 


A, pU2, drag coefficient 


Nomenclature of S. W. Churchill 


Ap projected area, m? 

Cy= aoe mean drag coefficient due to friction 
C= Teor , mean drag coefficient due to pressure 
C= ee total mean drag coefficient 

Fy drag force due to friction, N 

ya drag force due to pressure, N 


Bernoulli Theorems 


When the viscous effects are negligible compared with the inertial forces (i.e., large Reynolds 
number) there are a number of generalizations that can be made about the flow. These are 
described by the Bernoulli theorems. The fluid is assumed to be inviscid and have zero thermal 
conduction so that the flow is also barotropic (density a single-valued function of pressure). The 
first of the Bernoulli theorems is derived for flow that may be rotational. A special case is for 
motions relative to a rotating coordinate system where Coriolis forces arise. For irrotational 
flow, the Bernoulli theorem is a statement of the conservation kinetic energy, potential energy, 
and the expansion energy. A macroscopic energy balance can be made that includes the effects 
of viscous dissipation and the work done by the system. 


equations of motion for a Newtonian fluid is 
Equation: 


D 
or = pf—Vp+(A+p)V(Vev) + pV?v 


The assumptions of steady, inviscid flow simplify the equations to 
Equation: 


p(ve V)v = pf — Vp 


The assumptions of barotropic flow with conservative body forces allow, 
Equation: 


and by virtue of the identity 
Equation: 


(veV)v=V (Sv?) +wxv 


the equations of motion can be written 
Equation: 


VH=vxw 


If the body force is gravitational then 2 = pgz. 


Let H denote the function of which the gradient occurs on the left-hand side of this equation. 
VH is a vector normal to the surfaces of constant H. However, v x w is a vector 
perpendicular to both v and w so that these vectors are tangent to the surface. However, v and 
w are tangent to the streamlines and vortex lines respectively so that these must lie in a surface 
of constant H. It follows that H is constant along the streamlines and vortex lines. The surfaces 
of constant H which are crossed with this network of stream and vortex lines are known as 
Lamb surfaces and are illustrated in Fig. 6.2. 


Coriolis force. Suppose that the motion is steady relative to a steadily rotating axis with an 
angular velocity, w . Batchelor (1967) derives the equations of motion in this rotating frame and 
shows that we must now include the potential from the centrifugal force and the addition of a 
Coriolis force term. 

Equation: 


VH =v x (w+ 2w) 
2 
H=2+6(p)+iv-26 


Irrotational flow. If the flow is also irrotational, then w = 0 and hence the energy function 


Equation: 


H=2+6(p)+ 5v’ 


is constant everywhere. 


Ideal gas. For an ideal gas we have 
Equation: 


p= (cp — cv) pT 
and if the heat capacities are assumed constant an isentropic change of state results in the 


following expression for H. 
Equation: 


1.2 
Fiigeal gas — Q+ CpT oF 2U 
The gas is hotter at places on a streamline where the speed is smaller or has a lower potential 


energy. This may represent the heating of air at a stagnation point, cooling of ascending air, or 
heating of descending air. 


The transformation from the function of pressure ®(p) to heat capacity and temperature in the 
above equation for the isoentropic expansion of an ideal gas with constant heat capacity is 
proven as follows. The total differential of entropy in terms of pressure and temperature is as 


follows. 
Os Os 
dS = (e) ine ob (sr) a 


Equation: 
In an isoentropic expansion, the change of entropy is zero and this gives us a relation between 
the differential of pressure and differential of temperature. 


Equation: 
Os Os 
— =-—|——]dT, f = 
(5), (sr) 4 , for dS=0 


The coefficient on the left-hand side can be determined from the Maxwell relations and the 
ideal gas EOS. 
Equation: 


Maxwell relation 


—-~ 
TE 
n—"" 
mi 
| 
| 
— 
aS 
wn” 


1 . 
; OT? ideal gas 


The coefficient of the right hand side can be determined from the definition of the heat capacity 
at constant pressure. 


Equation: 
C5. = (33) , reversible process 
Pp 
= Os 
= 2 (sr). 
(9%) =2 
‘\ OT /p T 


Substituting these relations into the equation gives us the following. 


Equation: 
dp _ 
ps => Cp aT 
a T 
J i J cp aT 
= ¢)T, if cp is constant 


Q.E.D. 


Solution of the Partial Differential Equations 
Topics Covered in the Chapter 


¢ Classes of partial differential equations 

e Systems described by the Poisson and Laplace equation 
e Systems described by the diffusion equation 

¢ Greens function, convolution, and superposition 

¢ Green's function for the diffusion equation 

e Similarity transformation 

¢ Complex potential for irrotational flow 

¢ Solution of hyperbolic systems 


Classes of partial differential equations 


The partial differential equations that arise in transport phenomena are usually the first 
order conservation equations or second order PDEs that are classified as elliptic, 
parabolic, and hyperbolic. A system of first order conservation equations is sometimes 
combined as a second order hyperbolic PDE. The student is encouraged to read R. 
Courant, Methods of Mathematical Physics, Volume II Partial Differential Equations, 
1962 for a complete discussion. 


System of conservation laws. Denote the set of dependent variables (e.g., velocity, 
density, pressure, entropy, phase saturation, concentration) with the variable u and the 
set of independent variables as t and x, where x denotes the spatial coordinates. In the 
absence of body forces, viscosity, thermal conduction, diffusion, and dispersion, the 
conservation laws (accumulation plus divergence of the flux and gradient of a scalar) 
are of the form 


Equation: 
dg(u) _ Ogu) | Of(u) _ 9 af a 
pp (a ge Ge Pg Be 
where 
Of — O(g1, 925° : -gn) “a (fi, fo, "act +, fn) : s 
an | sens ~ sa) Gia ae Jacobian matrix 


This is a system of first order quasilinear hyperbolic PDEs. They can be solved by the 
method of characteristics. These equations arise when transport of material or energy 
occurs as a result of convection without diffusion. 


The derivation of the equations of motion and energy using convective coordinates 
(Reynolds transport theorem) resulted in equations that did not have the accumulation 
and convective terms in the form of the conservation laws. However, by derivation of 


the equations with fixed coordinates (as in Bird, Stewart, and Lightfoot) or by 
application of the continuity equation, the momentum and energy equations can be 
transformed so that the accumulation and convective terms are of the form of 
conservation laws. Viscosity and thermal conductivity introduce second derivative 
terms that make the system non-conservative. This transformation is illustrated by the 
following relations. 


Equation: 

i pVev = #24V =0 inui 

pi tepVev = 3 +Ve(pv)=0, continuity equation 
pe = p= + pve VF, for F=v, S, or E 
OApF) __ OF Op 
ae — Pa tha 

Ve(pFv) = pFVeviveV(pF)=pFVeviFveVp+pve VF 
pee = OF) Fe + Ve(pFv)—pFVev—FveVp 
= oF) + Ve (oF v) —F| oe +pVev+ve Vol 
O(pF 

p2e = 207) + Ve (p Fv), Q.E.D. 


Assignment 7.1 


a. For the case of inviscid, nonconducting fluid, in the absence of body forces, derive 
the steps to express the continuity equation, equations of motion, and energy 
equations as conservation law equations for mass, momentum, the sum of kinetic 
plus internal energy, and entropy. 

b. For the case of isentropic, compressible flow, express continuity equation and 
equations of motions in terms of pressure and velocity. Transform it to a second 
order hyperbolic equation in the case of small perturbations. 


Second order PDE. The classification of second order PDEs as elliptic, parabolic, and 
hyperbolic arise from a transformation of the independent variables. The classification 
apply to quasilinear (i.e., linear in the highest order derivatives) but we will only 
discuss linear equations with constant coefficients here. Numerical solutions are needed 
for quasilinear systems. Again let u denote the dependent variables and t, x, y, z as the 
independent variables. Examples of the different classes of equations are 

Equation: 


0= ou o Oru t p= V7u-+p, elliptic equation 


Ou _ Ou | Ou , Ou 


= a op or eS V-u+p, parabolic equation 


2 2 2 2 . . 
we = oe oo oe t+ p= V’u+p, hyperbolic equation 


The p term represents sources. When the cgs system of units is used in electrostatics 
and p is the charge density, the source is expressed as 47 p . The factor 47 is absent 
with the mks or SI system of units. The parabolic PDEs are sometimes called the 
diffusion equation or heat equation. In the limit of steady-state conditions, the parabolic 
equations reduce to elliptic equations. The hyperbolic PDEs are sometimes called the 
wave equation. A pair of first order conservation equations can be transformed into a 
second order hyperbolic equation. 


Convective-diffusion equation. The above equations represented convection without 
diffusion or diffusion without convection. When both the first and second spatial 
derivatives are present, the equation is called the convection-diffusion equation. 
Equation: 


du du 1 02u 
Ot Ox Noe Ox2 


Usually a dimensionless group such as the Reynolds number, or Peclet number appears 
as a factor to quantify the relative contribution of convection and diffusion. 


Systems described by the Poisson and Laplace equation 


We saw earlier that an irrotational vector field can be expressed as the gradient of a 
scalar and if in addition the vector field is solenoidal, then the scalar potential is the 
solution of the Laplace equation. 


Equation: 
v=-V69, irrotational flow 
Vev=O0=-V’¢ 
Vev=0=~—V*¢, incompressible, irrotational flow 


Also, if the velocity field is solenoidal then the velocity can be expressed as the curl of 
the vector potential and the Laplacian of the vector potential is equal to the negative of 
the vorticity. If the flow is irrotational, then the vorticity is zero and the vector potential 
is a solution of the Laplace equation. 


Equation: 


v=VxA, incompressible flow 

Vxv=w=-V’A, for VeA=0 

Vem = —w, in two dimensions 

Vxv=0=—-V*A, irrotational flow and VeA=0 

V-w = 0, for two dimensional, irrotational, incompressible flow 


Other systems, which are solution of the Laplace equation, are steady state heat 
conduction in a homogenous medium without sources and in electrostatics and static 
magnetic fields. Also, the flow of a single-phase, incompressible fluid in a homogenous 
porous media has a pressure field that is a solution of the Laplace equation. 


Systems described by the diffusion equation 


Diffusion phenomena occur with viscous flow, thermal conduction, and molecular 
diffusion. Heat conduction and diffusion without convection are described by the 
diffusion equation. Convection is always present in fluid flow but its contribution to the 
momentum balance is neglected for creeping (low Reynolds number) flow or cases 
where the velocity is perpendicular to the velocity gradient. In this case 

Equation: 


oe ai = “p +vV?v, velocity perpendicular to velocity gradient 


as —vV’v, if f and Vp vanish 


Green's function, convolution, and superposition 


A property of linear PDEs is that if two functions are each a solution to a PDE, then the 
sum of the two functions is also a solution of the PDE. This property of superposition 
can be used to derive solutions for general boundary, initial conditions, or distribution 
of sources by the process of convolution with a Green's function. The student is 
encouraged to read P. M. Morse and H. Feshbach, Methods of Theoretical Physics, 
1953 for a discussion of Green's functions. 


The Green's function is used to find the solution of an inhomogeneous differential 
equation and/or boundary conditions from the solution of the differential equation that 
is homogeneous everywhere except at one point in the space of the independent 


variables. (The initial condition is considered as a subset of boundary conditions here.) 
When the point is on the boundary, the Green's function may be used to satisfy 
inhomogeneous boundary conditions; when it is out in space, it may be used to satisfy 
the inhomogeneous PDE. 


The concept of Green's solution is most easily illustrated for the solution to the Poisson 
equation for a distributed source p(x, y, z) throughout the volume. The Green's function 
is a solution to the homogeneous equation or the Laplace equation except at (0, Yo, Zo) 
where it is equal to the Dirac delta function. The Dirac delta function is zero 
everywhere except in the neighborhood of zero. It has the following property. 
Equation: 


iy. ” f(Qs(é— 2) dé = f (2) 


The Green's function for the Poisson equation in three dimensions is the solution of the 
following differential equation 
Equation: 


V’G = —6 (x — x5) = —5(@ — 24)5(y — yo) 5 (z— 2) 
G (x|x.) = 7 


4n|x—Xo| 


It is a solution of the Laplace equation except at x = x, where it has a singularity, i.e., 
it has a point source. The solution of the Poisson equation is determined by convolution. 
Equation: 


uo) = ff [ Cx.) 0%) dts dys dz. 


Suppose now that one has an elliptic problem in only two dimensions. One can either 
solve for the Green's function in two dimensions or just recognize that the Dirac delta 
function in two dimensions is just the convolution of the three-dimensional Dirac delta 
function with unity. 

Equation: 


S(t ~ 0) 5(y-vo) = [5 (2-25) 5(y~ Yo) 5 (2 ~ 2) dee 


Thus the two-dimensional Green's function can be found by convolution of the three 
dimensional Green's function with unity. 
Equation: 


G(x, y|Zo, Yo) ee G(x|x,) dz, 


— glu |(w — 20)” + (y- yo)"| 


This is a solution of the Laplace equation everywhere except at (9, yo) where there is a 
line source of unit strength. The solution of the Poisson equation in two dimensions can 
be determined by convolution. 

Equation: 


u(z, y) _ // G(z, y|Lo, Yo) P(Lo, Yo) dx_ dYo 


Assignment 7.2 Derivation of the Green's function 

Derive the Green's function for the Poisson equation in 1-D, 2-D, and 3-D by 
transforming the coordinate system to cylindrical polar or spherical polar coordinate 
system for the 2-D and 3-D cases, respectively. Compare the results derived by 
convolution. 


Method of Images 


Green's functions can also be determined for inhomogeneous boundary conditions (the 
boundary element method) but will not be discussed here. The Green's functions 
discussed above have an infinite domain. Homogeneous boundary conditions of the 
Dirichlet type (u = 0) or Neumann type (Ou/On = 0) along a plane(s) can be 
determined by the method of images. 


Suppose one wished to find the solution to the Poisson equation in the semi-infinite 
domain, y > 0 with the specification of either uw = 0 or Ou/On = 0 on the boundary, 
y = 0. Denote as u? (2, y, z) the solution to the Poisson equation for a distribution of 
sources in the semi-infinite domain y > 0. The solutions for the Dirichlet or Neumann 
boundary conditions at y = 0 are as follows. 

Equation: 


u(x, y, 2) = yu? (x, y, Z) er Ts (z,—-y,z), forw=Oaty=0 
u(z,y, z) =u" (x,y, z) + ul (z, —Y; z), for du/dy = Oaty=0 


The first function is an odd function of y and it vanishes at y = 0. The second is an 
even function of y and its normal derivative vanishes at y = 0. 


An example of the method of images to satisfy either the Dirichlet or Neumann 
boundary conditions is illustrated in the following figure. The black curve is the 
response to a line sink at x = 1.5. We desire to have either the function or the 
derivative at x = 0 to vanish. The red curve is a line sink at x = —1.5. The sum of the 
two functions is symmetric about xz = 0 and has zero derivative there. The difference is 
anti-symmetric about z = 0 and vanishes at x = 0. 


Method of Images 


Now suppose there is a second boundary that is parallel to the first, i.e. y = a that also 
has a Dirichlet or Neumann boundary condition. The domain of the Poisson equation is 
now 0 < y < a. Denote as u! the solution that satisfies the BC at y = 0. A solution 
that satisfies the Dirichlet or Neumann boundary conditions at y = a are as follows. 
Equation: 


u(a,y,z) =u' (x,y,z) —ul(x,2a—y,z), foru=Oaty=a 


u(a,y,z) =u! (x,y,z) +ul(x,2a—y,z), for du/dy=Oaty=a 


This solution satisfies the solution at y = a, but no longer satisfies the solution at y = 0 
. Denote this solution as u? and find the solution to satisfy the BC at y = 0. By 
continuing this operation, one obtains by induction a series solution that satisfies both 
boundary conditions. It may be more convenient to place the boundaries symmetric 
with respect to the axis in order to simplify the recursion formula. 


Assignment 7.3 

Calculate the solution for a unit line source at the origin of the z,y plane with zero flux 
boundary conditions at y= +1 and y = —1. Prepare a contour plot of the solution for 
0 < x < 5. What is the limiting solution for large x? Note: The boundary conditions 
are conditions on the derivative. Thus the solution is arbitrary by a constant. 


Existence and Uniqueness of the Solution to the Poisson Equation 


If the boundary conditions for Poisson equation are the Neumann boundary conditions, 
there are conditions for the existence to the solution and the solution is not unique. This 
is illustrated as follows. 

Equation: 


Vu=—p inV, neVu=fonS 


If Vuav =~ [ffoav 
GnVuds =— fff pdv 
Gras =— [ff pav 


This necessary condition for the existence of a solution is equivalent to the statement 
that the flux leaving the system must equal the sum of sources in the system. The 
solution to the Poisson equation with the Neumann boundary condition is arbitrary by a 
constant. If a constant is added to a solution, this new solution will still satisfy the 
Poisson equation and the Neumann boundary condition. 


Green's function for the diffusion equation 


We showed above how the solution to the Poisson equation with homogeneous 
boundary conditions could be obtained from the Green's function by convolution and 
method of images. Here we will obtain the Green's function for the diffusion equation 
for an infinite domain in one, two, or three dimensions. The Green's function is for the 
parabolic PDE 

Equation: 


where the parameter a? represents the ratio of the storage capacity and the conductivity 
of the system and p is a known distribution of sources in space and time. The infinite 
domain Green's function gy, (x, tao, t.) is a solution of the following equation 
Equation: 


2 99n(X, t|Xo,to) 
a — 


V7 9n(X, t]Xo, to 
Gn(Xst|Xo, te) ai 


5 (x — xo) 5(t — te) 


The source term is an impulse in the spatial and time variables. The form of the Green's 
function for the infinite domain, for n dimensions, is (Morse and Feshbach, 1953) 
Equation: 


cab. if ca ‘ ~(a? R?/4r) +0 
Gn(X, t|Xo, to) = gn (FR, 7) = a? (ss) @ x 
0, T=0 
where 
T=t-—t, 


R=|x—-x,| 


This Green's function satisfies an important integral property that is valid for all values 
of n: 
Equation: 


1 
[ontRiryav Fa T>O0 


This expression is an expression of the conservation of heat energy. At a time t, at Xo, a 
source of heat is introduced. The heat diffuses out through the medium, but in such a 
fashion that the total heat energy is unchanged. 


The properties of this Green's function can be more easily seen be expressing it in a 
standard form 
Equation: 


a)" -[R2/2(27/0*)] S0 
029, (R,t) = (Es) ° _ 
0 , T<O0 


The normalized function ag, for n = 3 represent the probability distribution of the 
location of a Brownian particle that was at x, at time t,. The cumulative probability is 
equal to unity. 


The same normalized function for n = 1, corresponds to the normal or Gaussian 
distribution with the standard deviation given by 
Equation: 


_ Vir 


a 


on 


Observe the Green's function in one, two, and three dimensions by executing greens.m 
and the function, greenf.m in the diffuse subdirectory of CENG501. You may wish to 
use the code as a template for future assignments. 


Step Response Function The infinite domain Green's function is the impulse response 
function in space and time. The response for a distribution of sources in space or as an 
arbitrary function of time can be determined by convolution. In particular the response 
to a constant source for 7 > 0 is the step response function. It has classical solutions in 
one and two dimensions. The unit step function or Heaviside function is the integral of 
the Dirac delta function. 

Equation: 


1. ES SO 
6(t'—t,)dt'’=S(t—t,) =< i 


The response function to a unit step in the source can be determined by integrating the 
Greens function or the impulse response function in time. 
Equation: 


t t 
J [Von — 02.4] de’ = — 6 (xe — ses) f(t ~ te) ae 


—co 


t a(_f ond’) 
m( f on dt’) = 0° SS = OK = XG) SEE) 
VU, — a? S = —§(x—x,) S(t — to) 
where 


t 
On =f gaat 


In one dimension, the step response function that has a unit flux at x = 0 is (R. V. 
Churchill, Operational Mathematics, 1958) (note: source is 26(x — 0) S(t — 0)) 


Equation: 


= t =e 
Ufiex=—1 — 9974/ —— ew? — a? |x| erfc eve t>0 
? 7a? 


/ 4t/a? 


For comparison, the function that has a value of unity at z = 0 (Dirichlet boundary 
condition) is 
Equation: 


In two dimensions, the unit step response function for a continuous line source is (H. S. 
Carslaw and J. C. Jaeger, Conduction of Heat in Solids, 1959) 
Equation: 


= — =e & Hi( oe), t>0 
R= @=a,); ys ¥o)- 
—Fi(-z) = f° —du 


x U 


= expint (zx), MATLAB function 


For large times this function can be expressed as 
Equation: 


2 a . 
parr = a? ip (s ) —%, for 4 > 100 


+ = 0.5772... 


In three dimensions, the unit step response function for a continuous point source is (H. 
S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, 1959) 
Equation: 


Note:The a? .factor has the units of time/Z?. If time is made dimensionless with 
respect to a”/R? and R with respect to R, , then the factor will disappear from the 
argument of the erfc. 


Assignment 7.4 

Plot the profiles of the response to a continuous source in 1, 2, and 3 dimensions using 
the MATLAB code contins.m and continf.m in the diffuse subdirectory. From the 
integral of the profiles as a function of time, determine the magnitude, spatial and time 
dependence of the source. Note: The exponential integral function, expint will give 
error messages for extreme values of the argument. It still computes the correct values 
of the function. 


Convective-Diffusion Equation 


The convective-diffusion equation in one dimension will be expressed in terms of 

velocity and dispersion, 

Equation: 
Ou du _ fv O7u 
a + Va, = Kan 
we w= Oy Se 0 
u(0,t)=1, t>0 


The independent variables can be transformed from (, t) to a spatial coordinate that 
translates with the velocity of the wave in the absence of dispersion, (y, t). 
Equation: 


y=axrx-—vt 


This transforms the equation to the diffusion equation in the transformed coordinates. 
Equation: 


To see this, we will transform the differentials from x to y. 
Equation: 

Oy _ 

a. —v 


Oy __ 
rae 


The total differentials expressed as a function of (z, t) or (y, ) are equal to each other. 
Equation: 


du = (3) dt + (32) de 
du = (4) t+ (%) dy 


The total differentials expressed either way are equal. The partial derivatives in t and z 
can be expressed in terms of partial derivatives in t and y by equating the total 
differentials with either dt or dx equal to zero and dividing by the non-zero differential. 
Equation: 
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LO ON OO 
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Soa Se 
cK 
a. 
2 


Q/|& 
ah 
NS 

cb 


Substitution into the original equation results in the transformed equation. This result 
could have been derived in fewer steps by using the chain rule but would not have been 
as enlightening. 


The boundary condition at z = 0 is now at changing values of y. We will seek an 
approximate solution that has the boundary condition u(y — —oo) = 1. A simple 
solution can be found for the following initial and boundary conditions. 
Equation: 


1, y<0 
eis y=0 

(0, y>0 
u(y > —oo,t) =1 
u(y > oo, t) = 0 


This system is a step with no dispersion at t = 0. Dispersion occurs for t > 0 as the 
wave propagates through the system. The solution can be found with a similarity 
transform, which we will discuss later. For now, the approximate solution is given as 


Equation: 
—vut 
u= } erfe( y ) = berfe( = i ) 
JAK t JAK t 


The boundary condition at z = 0 will be approximately satisfied after a small time 
unless the Peclet number is very small. 


Similarity transformation 


In some cases a partial differential equation and its boundary conditions (and initial 
condition) can be transformed to an ordinary differential equation with boundary 
conditions by combining two independent variables into a single independent variable. 
We will illustrate the approach here with the diffusion equation. It will be used later for 
hyperbolic PDEs and for the boundary layer problems. 


The method will be illustrated for the solution of the one-dimensional diffusion 
equation with the following initial and boundary conditions. The approach will follow 
that of the Hellums-Churchill method. 
Equation: 

& = K24, t>0,2>0 

u(x,0) = ure 

u (0,t) = usc 


The PDE, IC and BC are made dimensionless with respect to reference quantities. 
Equation: 


*x* _ Uc 
U ‘le 
eo CoS) he 
t "a 
ie 
Lo 
OuX __ Kt, | 0?u* 
ot x2 Ox*? 
“(eo 0) = -.0 
u*(0, t*) | seem | 1, = Uo = UBC — UIC 


There are three unspecified reference quantities and two dimensionless groups. The BC 
can be specified to equal unity. However, the system does not have a characteristic time 
or length scales to specify the dimensionless group in the PDE. This suggests that the 
system is over specified and the independent variables can be combined to specify the 
dimensionless group in the PDE to equal 1/4. 


Equation: 
Kt 1 x : : : 
o/—lf wy is dimensionless 
x2 | 4 ul VA4K t 


u(x,t) = u(n) 


The partial derivatives will now be expressed as a function of the derivatives of the 
transformed similarity variable. 


Equation: 
a 
Ox V4Kt 
OF ss ER 
at 2t 
Ou _ du On _ —n du 
Ot ~— dn Ot ~~ 2 dy 
ou 1 du 
Ox J4Kt an 
Ou __ 1 du 
Oz? AK ts dr? 


The PDE is now transformed into an ODE with two boundary conditions. 
Equation: 


d?u* du* __ 
ae + ONG, = 9 
u*(n = 0) =1 
u*(n > co) = 0 


Equation: 


v=C, eo? 

ut= OC, fe" dn+ Co 
u*(n=0)=1 > Ch=1 

u*(n + 00) =0=Cy fp” eT dn+1 


GQ == 
1 for edn 


foe” dn 
* — —Jo © = 
Oe aig, P= ee) 


a (et) = erfe( Z| 


Therefore, we have a solution in terms of the combined similarity variable that is a 
solution of the PDE, BC, and IC. 


Complex potential for irrotational flow 


Incompressible, irrotational flows in two dimensions can be easily solved in two 
dimensions by the process of conformal mapping in the complex plane. First we will 
review the kinematic conditions that lead to the PDE and boundary conditions. Because 
the flow is irrotational, the velocity is the gradient of a velocity potential. Because the 
flow is solenoidal, the velocity is also the curl of a vector potential. Because the flow is 
two dimensional, the vector potential has only one non-zero component that is 
identified as the stream function. The kinematic condition at solid boundaries is that the 
normal component of velocity is zero. No condition is placed on the tangential 
component of velocity at solid surfaces because the fluid must be inviscid in order to be 
irrotational. 

Equation: 


0g O 
v = Vo= (32, 3,0) = (Up, y, Uz) 
Vev = 0=V'¢ 
v= Vx A= (SS, es. 0) 


) -—d 
= (3, sui) 0) = (Vz, Vy, Vz) 


d¢ dg = Op —dp 
(=. 3,0) = (3, a 0) 
06 _ = Oy 
Or t: ~ oy 
06 _ = - Oy 
and Oy = Oz 


Both functions are a solution of the Laplace equation, i.e., they are harmonic and the 
last pair of equations corresponds to the Cauchy-Riemann conditions if ¢ and w are the 
real and imaginary conjugate parts of a complex function, w(z). 

Equation: 


wz) = o(z)+i¥(z) 
ZS EAST 


1/2 


= re“=r(cosO+isin@), r=|z|=(z?+y’) 6 =arctan (y/x) 


¢(z) = real[w(z)| 
W(z) imaginary|w(z)| 


The Cauchy-Riemann conditions are the necessary and sufficient condition for the 
derivative of a complex function to exist at a point z, , i.e., for it to be analytical. The 
necessary condition can be illustrated by equating the derivative of w(z) taken along 
the real and imaginary axis. 

Equation: 


aul) Re(w! (z)) + iIm(w’ (2) 
— lim setieh _ 8 get 


dt+i0 Oa 
=o bpt+idy 4 _ Ow - Ob 
— O+idy i” Oy 1 Oy 
Og _ Op 
Ox ~~ ~Oy 
ob _ _ ob 
and 57 = Dy? Q.E.D. 


Also, if the functions have second derivatives, the Cauchy-Riemann conditions imply 
that each function satisfies the Laplace equation. 


Equation: 
#6 | 0% 9 
Ox? Oy 
ay | Wb _ 4g 
Ox? Oy2 


The Cauchy-Riemann conditions also imply that the gradient of the velocity potential 
and the stream function are orthogonal. 
Equation: 


Ag Ob | Od Ob _ 


Wa yea) = Ox Ox | Oy Oy 


If the gradients are orthogonal then the equipotential lines and the streamlines are also 
orthogonal, with the exception of stagnation points where the velocity is zero. 


Since the derivative 
Equation: 


Os oe 
dz |dz—-0| 6z 


is independent of the direction of the differential 6z in the (a, y) -plane, we may 
imagine the limit to be taken with 6z remaining parallel to the x-axis (6z = 6x) giving 
Equation: 

dw Ob ,OY | 


= 1 = 92 —40 
dz Ox Ox if 


Now choosing z to be parallel to the y-axis (dz = idy), 
Equation: 
dw 1 0¢ , Ow 


dz i oy Oy 


These equations are a restatement that an analytical function has a derivative defined in 
the complex plane. Moreover, we see that the real part of w’ (z) is equal to v, and the 
imaginary part of w’ (z) is equal to —vy. If v is written for the magnitude of v and for 
the angle between the direction of v and the x-axis, the expression for dw/dz becomes 
Equation: 


dw. 20 


dz 
or 


= Uz —1Vy = VE 


i= real | @ | 


Vy = — imaginary | 4 | 


Flow Fields. The simplest flow field that we can imagine is just a constant translation, 
w = (U — iV )z where U and V are real constants. The components of the velocity 
vector can be determined from the differential. 

Equation: 


w(z) = (U -iW)z=(U -wW)(a+iy) =Ur2+Vy+i(-V2+Uy)=¢+in 
m@ = (U-iV) =v, —idy 

vz =U, vw=aV 

g=Uzr+Vy, w=-Vzxe+Uy 


Another simple function that is analytical with the exception at the origin is 
Equation: 


w(z) = Az™=Arre'r? 
= Ar”"cos nO6+iAr” sin nO 
= g+iy 
thus 
@ = Ar” cos no 
wy = Ar”sin nO 
dw 


ee VAg Si 


where A and n are real constants. The boundary condition at stationary solid surfaces 
for irrotational flow is that the normal component of velocity is zero or the surface 
coincides with a streamline. The expression above for the stream function is zero for all 
r when 6 = 0 and when 0 = z/n. Thus these equations describe the flow between 
these boundaries are illustrated below. 


as 


Fig. 6.5.1 (Batchelor, 1967) Irrotational flow in 
the region between two straight zero-flux 
boundaries intersecting at an angle 7/n. 


Earlier we discussed the Green's function solution of a line source in two dimensions. 
The same solution can be found in the complex domain. A function that is analytical 
everywhere except the singularity at z, is the function for a line source of strength m. 
Equation: 


w(z) = = In(z—2z,), line source 
dw. Fe aN te 
dz 2m (z-Z) VE hy 


This results can be generalized to multiple line sources or sinks by superposition of 
solutions. A special case is that of a source and sink of the same magnitude. 
Equation: 


w(z) = >0,5¢ In(z- 2), multiple line sources 


_ om Z—Zo a : 
w(z)= 3 In (= ae ) source — sink pair 


The above flow fields can be viewed with the MATLAB code corner.m, linesource.m, 
and multiple.m in the complex subdirectory. 


Assignment 7.5 

Line Source Solution For 2, at the origin, derive expressions for the flow potential, 
stream function, components of velocity, and magnitude of velocity for the solution to 
the line source in terms of r and @ . Plot the flow potentials and stream functions. 
Compute and plot the flow potentials and stream function for the superposition of 
multiple line sources corresponding to the zero flux boundary conditions at y = +1 and 
—1 of the earlier assignment. 


The circle theorem. (Batchelor, 1967) The following result, known as the circle theorem 
(Milne-Thompson, 1940) concerns the complex potential representing the motion of an 
inviscid fluid of infinite extent in the presence of a single internal boundary of circular 
form. Suppose first that in the absence of the circular cylinder the complex potential is 
Equation: 


w° = f(z) 


and that f(z) is free from singularities in the region |z| < a, where a is a real length. If 
now a Stationary circular cylinder of radius a and center at the origin bounds the fluid 
internally, the flow is modified to the following complex potential: 

Equation: 


w! = f(z) + F(a?/2) 


We show that the surface of the cylinder, |z| = a, is a streamline. 
Equation: 


Equation: 


w*(z)| ig = f@2D)ljjeat+ f (2/2) pe 
= f(Dlagcat FD, 
= f@Djgoat+ f @ ey 
- 2 Real] f(2)|j.-4 +40 
= Plzl=a +t Yiz=c 


A complex potential of this form thus has |z| = a as a streamline, ~ = 0; and it has the 
same singularities outside |z| = a as f(z), since if z lies outside |z| = a, a? /z lies in 
the region inside this circle where f(z) is known to be free from singularities. 
Consequently the additional term f. (a? / z) in the equation represents fully the 
modification to the complex potential due to the presence of the circular cylinder. It 
should be noted that the complex potentials considered, both in the absence and in the 
presence of the circular cylinder, refer to the flow relative to axis such that the cylinder 
is stationary. 


The simplest possible application of the circle theorem is to the case of a circular 
cylinder held fixed in a stream whose velocity at infinity is uniform with components 
(—U, —V). In the absence of the cylinder the complex potential is —(U — i V)z, it is 
singular at infinity and the circle theorem shows that, with the cylinder present, 
Equation: 


w(z) = —-(U-iV)z—(U+iV)a*/z 


The potentials and streamlines for the steady translation of an inviscid fluid past a 
circular cylinder can be viewed with the MATLAB code circle.m. 


Potentials and Streamlines for Flow Past Cylinder 


Conformal Transformation (Batchelor, 1967). We now have the complex potential flow 
solutions of several problems with fairly simple boundary conditions. These solutions 
are analytical functions whose real and imaginary parts satisfy the Laplace equation. 
They also have a streamline that coincides with the boundary to satisfy the condition of 
zero flux across the boundary. Conformal transformations can be used to obtain 
solutions for boundaries that are transformed to different shapes. Suppose we have an 
analytical function w(z) in the z = x + iy plane. This solution can be transformed to 
the ¢ = € + in plane as another analytical function provided that the relation between 
these two planes, ¢ = F(z) is an analytical function. This mapping is a connection 
between the shape of a curve in the z - plane and the shape of the curve traced out by 
the corresponding set of points in the ¢ — plane. The solution in the ¢ — plane is 
analytical, i.e., its derivative defined, because the mapping, ¢ = F'(z) is an analytical 
function. The inverse transformation is also analytical. 

Equation: 


dw(s) _ du(z) dz _ & 
ds ~—s dz ds ~~ 4F(z) 


dw(z) _ dw(s) dF(z) 
dz — dg dz 


w(C¢) is thus the complex potential of an irrotational flow in a certain region of the ¢- 
plane, and the flow in the z-plane is said to been 'transformed' into flow in the ¢-plane. 
The family of equipotential lines and streamlines in the z-plane given by 

g(x, y) = const. and w(x, y) = const. transform into families of curves in the ¢-plane 
on which ¢ and ~# are constant and which are equipotential lines and streamlines in the 
¢-plane. The two families are orthogonal in their respective plane, except at singular 
points of the transformation. The velocity components at a point of the flow in the ¢- 
plane are given by 

Equation: 


_; Se ae _; yee 
VE Wy = de = ae de = (Uz Wy de 


This shows that the magnitude of the velocity is changed, in the transformation from the 
z-plane to the ¢-plane, by the reciprocal of the factor by which linear dimensions of 
small figures are changed. Thus the kinetic energy of the fluid contained within a closed 
curve in the z-plane is equal to the kinetic energy of the corresponding flow in the 
region enclosed by the corresponding in the ¢-plane. 


outside of an ellipse in the z-plane into the region outside a circle in the ¢-plane is given 
by 
Equation: 

2 


ae es 


c= rz 5 (2? An?)r? 


where J is a real constant so that 


Equation: 
r? r? 
ao We ee ae A eek / | OF ceria 
I<| Is| 


This converts a circle of radius c with center at the origin in the - plane into the ellipse 


Equation: 


in the z-plane, where 
Equation: 


eer 


If the ellipse is mapped into a circle in the ¢-plane, it is convenient to use polar 
coordinates (r, @), especially since the boundary corresponds to a constant radius. The 
radius that maps to the elliptical boundary is (ellipse.m in the complex directory) 


Equation: 
i= : lo oe 
eo 2 : a—b 


The transformation from the polar coordinates to the z - plane is defined by 
Equation: 


z= 2X cosh w 
where 
w=r+id 


Confocal E llipses and Conjugate Hyperbolae, a = 1 b= 0.2 


Transformation of cylindrical polar 
coordinates into orthogonal, elliptical 
coordinates 


The polar coordinates (r, @), transform to an orthogonal set of curves which are 
confocal ellipses and conjugate hyperbolae. 


This transformation can be substituted into the complex potential expression for the 
flow of a fluid past a circular cylinder. 
Equation: 


1 
w= —5 (2 + b)[(U —iV)e*"" + (U + iV)e”?*] 
It is convenient to write - for the angle which the direction of motion of the flow at 


infinity makes with the x - axis so that 
Equation: 


U+iV = (U2 +V?)"e-2 


The complex potential now becomes 


Equation: 


w= —(U2+?)?(a +0) cosh (w — r, + ia) 


The corresponding velocity potential and stream function are 
Equation: 


a= —(U? +.V)"?(a +0) cosh (r — r,) cos (0+ a) 


p= —(U? +.V?)"?(a +0) sinh (r — r,) sin (8+ a) 


The velocity potentials and streamlines are illustrated below for flow past an elliptical 
cylinder (fellipse.m in the complex directory). Note the stagnation streamlines on either 
side of the body. These two stagnation points are regions of maximum pressure and 
result in a torque on the body. Which way will it rotate? 


Potentials and Streamlines for Flow Past Ellipse 


Flow past an ellipse of an inviscid fluid that is in 
steady translation at infinity. 


Pressure distribution. When an object is in a flow field, one may wish to determine the 
force exerted by the fluid on the object, or the ‘drag’ on the object. Since the flow field 
discussed here has assumed an inviscid fluid, it is not possible to determine the viscous 
drag or skin friction directly from the flow field. It is possible to determine the 'form 
drag' from the normal stress or pressure distribution around the object. However, one 
must be critical to determine if the calculated flow field is physically realistic or if some 
important phenomena such as boundary layer separation may occur but is not allowed 
in the complex potential solution. 


The Bernoulli theorems give the relation between the magnitude of velocity and 
pressure. We have assumed irrotational, incompressible flow. If in addition we assume 
the body force can be neglected, then the quantity, H, must be constant along a 
streamline. 

Equation: 


H= a =f iu = constant 


p= + pv +constant, since p isalso constant 


The pressure relative to some datum can be determined by the square of the magnitude 
of velocity. This is easily calculated from the complex potential. 
Equation: 


There are some theorems that facilitate the integration of pressure around bodies in the 
complex plane, but they will not be discussed here. The pressure and tangential velocity 
profiles for the inviscid flow around an object are needed for calculation of the viscous 
flow in the boundary layer between the solid boundary and the inviscid outer flow. 


Assignment 7.6 

Pressure profiles Calculate the pressure field or the square of the velocity field for the 
flow in or around a corner and the flow past a circular cylinder. Look at the expression 
for the corner flow. Under what conditions is there a flow singularity? Show the 
pressure or velocity squared pseudo-color for wall angles of 1/2, 7, 37/2, and 27 . 
Which cases are physically realistic and what do you think happens in the unrealistic 
cases? What is special about the pressure profile around the circular cylinder? What 


value of form drag will it predict? Is it realistic and if not, why not? Add the following 
code to the code for corner flow and flow around a circular cylinder. pause 


% Calculate pressure distribution from square of velocity 
(your code here to calculate pressure field) 

pcolor(x,y,Pp) 

colormap(hot ) 

shading flat 

axis image 


Solution of hyperbolic systems 


The conservation equations for material, momentum, and energy reduce to first order 
PDE in the absence of diffusivity, dispersion, viscosity, and heat conduction. In thin 
films, viscosity may be a dominant effect in the velocity profile normal to the surface 
but the continuity equation integrated over the film thickness will have only first order 
spatial derivatives unless the effects of interfacial curvature become important. In one 
dimension, the system of first order partial differential equations can be calculated by 
the method of characteristics [A. Jeffrey (1976), H.-K. Rhee, R. Aris, and N. R. 
Amundson (1986, 1987)]. Here we will only consider the case of a single dependent 
variable with constant initial and boundary conditions. Denote the dependent variable as 
S and the independent variables as x and t. The differential equation with its initial and 
boundary conditions are as follows. 

Equation: 


8s 4 95) 9, t>0,2>0 
S (x, 0) = Sic 
f (0,t) = fac 


The dependent variable can be normalized such that the initial condition is equal to zero 
and the boundary condition is equal to unity. Thus, henceforth it is assumed such a 
transformation has been made. The dependent variable will be called 'saturation' and the 
flux called ‘fractional flow' to use the nomenclature for multiphase flow in porous 
media. However, the dependent variable could be film thickness in film drainage or 
height of a free surface as in water waves. The PDE, IC, and BC are rewritten as 
follows. 

Equation: 


9s af(S) as _ 
+ “dS Ox = 0, 


S{2,0) — 0 
f(0,t) =1 


t>0, z«>0 


The differential, df /dS is easily calculated since there is only one independent 
saturation. If there were three or more phases this differential would be a Jacobian 
matrix. The locus of constant saturation will be sought by taking the total differential of 


ole): 
Equation: 
Os OS 
dS = ——dt+ —dz=0 
Ot " 0 
Equation: 
dx = 0S/dt 
care ~  AS/da 
d 
5(5) 
= Us 


This equation expresses the velocity that a particular value of saturation propagates 
through the system, i.e., the saturation velocity, vg is equal to the slope of the fractional 


flow curve. It is also the slope of a trajectory of constant saturation (i.e., dS = 0) in the 
(x,t) space. Since we are assuming constant initial and boundary conditions, changes 
in saturation originate at (2, t) = (0,0). From there the changes in saturation, called 
waves, propagate in trajectories of constant saturation. We assume that df/dS isa 
function of saturation and independent of time or distance. This assumption will result 
in the trajectories from the origin being straight lines if the initial and boundary 
conditions are constant. The trajectories can easily be calculated from the equation of a 
straight line. 

Equation: 


df(S) 
dS 


“(S)= t 


Wave: A composition (or saturation) change that propagates through the system. 


X X 


Spreading wave: A wave in which neighboring composition (or saturation) values 
become more distant upon propagation. 


Equation: 
dx ) 2 dx ) 
dt Sa dt Sp 


tt to>t1 


X X 


Indifferent waves: A wave in which neighboring composition (or saturation) values 
maintain the same relative position upon propagation. 
Equation: 


(=) -($) 
dt oan dt Sb 


a tt 


X X 


Step Wave: An indifferent wave in which the compositions change discontinuously. 


>t 


X X 


Self Sharpening Waves: A wave in which neighboring compositions (saturations) 
become closer together upon propagation. 


Equation: 
dz . dz 
dt Sa dt Sp 


Shock Wave: A wave of composition (saturation) discontinuity that results from a self 
sharpening wave. 


x 


Rule: Waves originating from the same point (e.g., constant initial and boundary 
conditions) must have nondecreasing velocities in the direction of flow. This is another 
way of saying that when several waves originate at the same time, the slower waves can 
not be ahead of the faster waves. If slower waves from compositions close to the initial 
conditions originate ahead of faster waves, a shock will form as the faster waves 


overtake the slower waves. This is equivalent to the statement that a sharpening wave 
can not originate from a point; it will immediately form a shock. 


Spreading Wave 
Wave does 
not exist 


= = 
— 


Equation: 
x dx df 
a Co a ae 


Mass Balance Across Shock 


We saw that sharpening wave must result in a shock but that does not tell us the velocity 
of a shock nor the composition (saturation) change across the shock. To determine these 
we must consider a mass balance across a shock. This is sometimes called an integral 
mass balance as opposed to the differential mass balance derived earlier for continuous 
composition (saturation) changes. 


Equation: 


Equation: 


Equation: 


Equation: 


$1 


p AAx(S2 = Si) 


Accumulation : 
u AAt( fo — fi) 


input — output : 


y AAz AS =uAAtAf 


aD, 


Af/AS is the cord slope of the f versus S curve between S; and $3. 


The conservation equation for the shock shows the velocity to be equal to the cord slope 
between S$; and S> but does not in itself determine S; and S». To determine S; and So, 
we must apply the rule that the waves must have non-decreasing velocity in the 
direction of flow. The following figure is a solution that is not admissible. This 
solution is not admissible because the velocity of the saturation values (slope) between 
the IC’ and $1 are less than that of the shock and the velocity of the shock (cord slope) 
is less than that of the saturation values immediately behind the shock. 


Not Admissible 
BC 


shock oe 
f \ 
iC S1 
o> 


This solution is admissible in that the velocity in nondecreasing in going from the BC 
to the IC’. However, it in not unique. Several values of $2 will give admissible 
solutions. Suppose that the value shown here is a solution. Also suppose that dispersion 
across the shock causes the presence of other values of S between $1 and $2. There are 
some values of S that will have a velocity (slope) greater than that of the shock shown 
here. These values of S will overtake $2 and the shock will go the these values of S. 
This will continue until there is no value of S that has a velocity greater than that of the 
shock to that point. At this point the velocity of the saturation value and that of the 
shock are equal. On the graphic construction , the cord will be tangent to the curve at 
this point. This is the unique solution in the presence of a small amount of dispersion. 


admissible but not unique 
S2 
BC 


IC 


S 


Composition (Saturation) Profile The composition (saturation) profile is the 
composition distribution existing in the system at a given time. 


X 


Composition or Flux History: The composition or flux appearing at a given point in 
the system, e.g.,z = 1. 


Summary of Equations 


The dimensionless velocity that a saturation value propagates is given by the following 


equation. 
dx df 
Gee 
dS=0 


Equation: 

With uniform initial and boundary conditions, the origin of all changes in saturation is 
at z = 0 andt = 0. If f(S) depends only on S and not on z or t, then the trajectories 
of constant saturation are straight line determined by integration of the above equation 
from the origin. 


Equation: 
Ms d 
«($) = ()as-ot = ay (S)t 
«(AS) = (FZ) ast = gt 


These equations give the trajectory for a given value of S or for the shock. By 
evaluating these equations for a given value of time these equations give the saturation 
profile. 


The saturation history can be determined by solving the equations for ¢ with a specified 
value of z,e.g.%2 = 1. 


Equation: 
t(AS) = =, e=1 
AS 
— zx = 
(5) = aq ee 


The breakthrough time, t 37, is the time at which the fastest wave reaches x = 1.0. 
The flux history (fractional flow history) can be determined by calculating the fractional 
flow that corresponds to the saturation history. 


Summary of Diagrams 


The relationship between the diagrams can be illustrated in a diagram for the 
trajectories. The profile is a plot of the saturation at t = t,. The history at x = 1.0 is 
the saturation or fractional flow at x = 1. In this illustration, the shock wave is the 
fastest wave. Ahead of the shock is a region of constant state that is the same as the 
initial conditions. 


« t=to B.C. 
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Laminar Flows with Dependence on One Dimension 
¢ Couette flow 


o Planar Couette flow 
© Cylindrical Couette flow 
o Planer rotational Couette flow 


e Hele-Shaw flow 
Poiseuille flow 


o Friction factor and Reynolds number 
o Non-Newtonian fluids 


e Steady film flow down inclined plane 
e Unsteady viscous flow 


© Suddenly accelerated plate 
o Developing Couette flow 


Reading Assignment 
Chapter 2 of BSL, Transport Phenomena 


One-dimensional (1-D) flow fields are flow fields that vary in only one spatial dimension in 
Cartesian coordinates. This excludes turbulent flows because it cannot be one-dimensional. 
Acoustic waves are an example of 1-D compressible flow. We will concern ourselves here 
with incompressible 1-D flow fields that result from axial or planar symmetry. Cartesian, 1-D 
incompressible flows do not have a velocity component (other than possibly a uniform 
translation) in the direction of the spatial dependence because of the condition of zero 
divergence. Thus the nonlinear convective derivative disappears from the equations of motion 
in Cartesian coordinates. They may not disappear with curvilinear coordinates. 

Equation: 


v = v(z3) 
Ov 

Vev=-0 => Des =0 
v3 (x3 0) 0 > U3 0 

Ou; 
ve VW = 1105; = 035- = 0 
Ov; 0; . 
pa =—Vp+pf—-VzeT=—Vp+ pitugr, j=1,2 


We can demonstrate that this relation may not apply in curvilinear coordinates by considering 
an example with cylindrical polar coordinates. Suppose that the only nonzero component of 
velocity is in the @ direction and the only spatial dependence is on the r coordinate. The radial 
component of the convective derivative is non-zero due to centrifugal forces. 


Equation: 


v = (0, ve(r), 0] 


[ve Vv]. = a2 


The flows can be classified as either forced flow resulting from the gradient of the pressure or 
the potential of the body force or induced flow resulting from motion of one of the bounding 
surfaces. 


Some flow fields that result in 1-D flow are listed below and illustrated in the following 
figure (Churchill, 1988) 


. Forced flow through a round tube 

. Forced flow between parallel plates 

. Forced flow through the annulus between concentric round tubes of different diameters 
. Gravitational flow of a liquid film down an inclined or vertical plane 

. Gravitational flow of a liquid film down the inner or outer surface of a round vertical 


tube 


. Gravitational flow of a liquid through an inclined half-full round tube 
. Flow induced by the movement of one of a pair of parallel planes 
. Flow induced in a concentric annulus between round tubes by the axial movement of 


either the outer or the inner tube 


. Flow induced in a concentric annulus between round tubes by the axial rotation of either 


the outer or the inner tube 


. Flow induced in the cylindrical layer of fluid between a rotating circular disk and a 


parallel plane 


. Flow induced by the rotation of a central circular cylinder whose axis is perpendicular to 


parallel circular disks enclosing a thin cylindrical layer of fluid 


. Combined forced and induced flow between parallel plates 
. Combined forced and induced longitudinal flow in the annulus between concentric 


round tubes 


. Combined forced and rotationally induced flow in the annulus between concentric round 


tubes 


4 One-Dimensional Laminar Flows 
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Geometry and conditions that produce 
one-dimensional velocity fields 
(Churchill, 1988) 


Couette Flow 


The flows when the fluid between two parallel surfaces are induced to flow by the motion of 
one surface relative to the other is called Couette flow. This is the generic shear flow that is 
used to illustrate Newton's law of viscosity. Pressure and body forces balance each other and 
at steady state the equation of motion simplify to the divergence of the viscous stress tensor 
or the Laplacian of velocity in the case of a Newtonian fluid. 


Planar Couette flow. (case 7). 
Equation: 


d. d?v; 
=p =0, j=1,2 
dz3 


dx3 


The coordinates system can be defined so that v = 0 at x3 = 0 and the 7 component of 
velocity is non-zero at 73 = L. 


Equation: 
vj =0, 23=0 
vwj=U;, #23=L, jg=1,2 
The velocity field is 
Equation: 
U; L3 


The shear stress can be determined from Newton's law of viscosity. 
Equation: 

du j U j 
L L ? 


Fe 1,2 


Cylindrical Couette flow. The above example was the translational movement of two planes 
relative to each other. Couette flow is also possible in the annular gap between two concentric 
cylindrical surfaces (cases 8 and 9) if secondary flows do not occur due to centrifugal forces. 
We use cylindrical polar coordinates rather than Cartesian and assume vanishing Reynolds 
number. The only independent variable is the radius. 

Equation: 


: L(r vr) =0 
[Ver], = 12 (r Trr) — *£, may not vanish if Reynolds number is high 


[Vez],= 43 (r?7,9) + “— =0, may not vanish if Reynolds number is high 


Tr 


Newtonian fluid 


The stress profile can be calculated by integration. 
Equation: 


179 = 7°76). = TT | 


T=11 T=T2 


TT rz = TT ¢2| = TT rz| 


TT, i i) 


The boundary conditions on velocity depend on whether the cylindrical surfaces move 
relative to each other as a result of rotation, axial translation, or both. 
Equation: 


v= 0. = 


Vv =(0,06.U),. PST 


The velocity field for cylindrical Couette flow of a Newtonian fluid is . 
Equation: 


a Us r "1 
a wine (t 2 | 


Uz 


Ue Youtea fray. log (r/r1) 


Planer rotational Couette flow. The parallel plate viscometer has the configuration shown in 
case 10. The system is not strictly 1-D because the velocity of one of the surfaces is a 
function of radius. Also, there is a centrifugal force present near the rotating surface but is 
absent at the stationary surface. However, if the Reynolds number is small enough that 
secondary flows do not occur, then the velocity at a given value of the radius may be 
approximated as a function of only the z distance in the gap. The differential equations at 
zero Reynolds number are as follows. 

Equation: 


[Ver], = 5* =0 


Oz 
07 v9 
Oz? 


=0, Newtonian fluid 


Suppose the bottom surface is stationary and the top surface is rotating. Then the boundary 
conditions are as follows. 
Equation: 


veg = 0, z—0 
va=2rrQ, z=L 


The stress and velocity profiles are as follows. 
Equation: 


Tz = To2(r) = —2nr Ny (r)/L 
vg =2arrQz/L, Newtonian fluid 


The stress is a function of the radius and if the fluid is non-Newtonian, the viscosity may be 
changing with radial position. 


Plane-Poiseuille and Hele-Shaw flow 


Forced flow between two stationary, parallel plates, case 2, is called plane-Poiseuille flow or 
if the flow depends on two spatial variables in the plane, it is called Hele-Shaw flow. The 
flow is forced by a specified flow rate or a specified pressure or gravity potential gradient. 
The pressure and gravitational potential can be combined into a single variable, P. 
Equation: 


—Vp+ pf =-VP 
where 
P=p+pgh 


The product gh is the gravitational potential, where g is the acceleration of gravity and h is 
distance upward relative to some datum. The pressure, p, is also relative to a datum, which 
may be the datum for h. 


The primary spatial dependence is in the direction normal to the plane of the plates. If there is 
no dependence on one spatial direction, then the flow is truly one-dimensional. However, if 
the velocity and pressure gradients have components in two directions in the plane of the 
plates, the flow is not strictly 1-D and nonlinear, inertial terms will be present in the equations 
of motion. The significance of these terms is quantified by the Reynolds number. If the flow 
is steady, and the Reynolds number negligible, the equations of motion are as follows. 
Equation: 


~~ 62; 0x3? j=1,2 


2a). 
0= oe + re - j=1,2, Newtonian fluid 


Let h be the spacing between the plates and the velocity is zero at each surface. 
Equation: 


v; = 0, r3=0,h p= eZ 


The velocity profile for a Newtonian fluid in plane-Poiseuille flow is 
Equation: 


h? OP [;23\2 23 
oy che pe ee | ee eg _ < < 
v; 5 a | (2) | jg=1,2, VA ach 


The average velocity over the thickness of the plate can be determined by integrating the 
profile. 
Equation: 


h? @P 


ee ee ee 
19 One 


This equation for the average velocity can be written as a vector equation if it is recognized 
that the vectors have components only in the (1, 2) directions. 
Equation: 


¥=-——VP, ¥=V¥(x1,22), VP=VP(q1, 22) 


If the flow is incompressible, the divergence of velocity is zero and the potential, P, is a 
solution of the Laplace equation except where sources are present. If the strength of the 
sources or the flux at boundaries are known, the potential, P, can be determined from the 
methods for the solution of the Laplace equation. 


We now have the result that the average velocity vector is proportional to a potential gradient. 
Thus the average velocity field in a Hele-Shaw flow is irrotational. If the fluid is 
incompressible, the average velocity field is also solenoidal can can be expressed as the curl 
of a vector potential or the stream function. The average velocity field of Hele-Shaw flow is 
an physical analog for the irrotational, solenoidal, 2-D flow described by the complex 
potential. It is also a physical analog for 2-D flow of incompressible fluids through porous 


media by Darcy's law and was used for that purpose before numerical reservoir simulators 
were developed. 


Poiseuille Flow 


Poiseuille law describes laminar flow of a Newtonian fluid in a round tube (case 1). We will 
derive Poiseuille law for a Newtonian fluid and leave the flow of a power-law fluid as an 
assignment. The equation of motion for the steady, developed (from end effects) flow of a 
fluid in a round tube of uniform radius is as follows. 


Equation: 
— _ OP 
_ Or 
OP 1 0 
= — Fe ewe (Tre) O<r<R 


0= 2 aE pie (2), Newtonian fluid 


The boundary conditions are symmetry at r = 0 and no slip at r = R. 
Equation: 


| =, Ov; =) 
Trz|\p—9 — LO; 0 


T= 


(=), Pah 


From the radial component of the equations of motion, P does not depend on radial position. 
Since the flow is steady and fully developed, the gradient of P is a constant. The z 
component of the equations of motion can be integrated once to derive the stress profile and 
wall shear stress. 


Equation: 
Tie (-2)5, 0<r<R 
re= (-$)4 


If the fluid is Newtonian, the equation of motion can be integrated once more to obtain the 
velocity profile and maximum velocity. 


Equation: 
R aP r " 
vy, = # (-8P)\1-(4 O<r<R 
z~ “Au ( az ) (zz) : , Newtonian fluid 
R? (-2) 
Uz, max Au Oz 


The volumetric rate of flow through the pipe can be determined by integration of the velocity 
profile across the cross-section of the pipe, i.e.,0 << r< Rand0 <6 < 2z. 
Equation: 


R* oP 
Q= 2 (-=). Newtonian fluid 


This relation is the Hagen-Poiseuille law. If the flow rate is specified, then the potential 
gradient can be expressed as a function of the flow rate and substituted into the above 
expressions. 


The average velocity or volumetric flux can be determined by dividing the volumetric rate by 
the cross-sectional area. 
Equation: 


-2). Newtonian fluid 


Before one begins to believe that the Hagen-Poiseuille law is a "law" that applies under all 
conditions, the following is a list of assumptions are implicit in this relation (BSL, 1960). 


a. The flow is laminar-NRe less than about 2100. 

b. The density is constant ("incompressible flow"). 

c. The flow is independent of time ("steady state"). 

d. The fluid is Newtonian. 

e. End effects are neglected-actually an "entrance length" (beyond the tube entrance) on the 
order of Le = 0.035D NRe is required for build-up to the parabolic profiles; if the 
section of pipe of interest includes the entrance region, a correction must be applied. The 
fractional correction introduced in either P or Q never exceeds Le/L if L > Le. 

f. The fluid behaves as a continuum-this assumption is valid except for very dilute gases or 


very narrow capillary tubes, in which the molecular mean free path is comparable to the 
tube diameter ("slip flow" regime) or much greater than the tube diameter ("Knudsen 
flow" or "free molecule flow" regime). 

g. There is no slip at the wall-this is an excellent assumption for pure fluids under the 
conditions assumed in ( f ). 


Friction factor and Reynolds number. Because pressure drop in pipes is commonly used in 
process design, correlation expressed as friction factor versus Reynolds number are available 
for laminar and turbulent flow. The Hagen-Poiseuille law describes the laminar flow portion 
of the correlation. The correlations in the literature differ when they use different definitions 
for the friction factor. Correlations are usually are usually expressed in terms of the Fanning 
friction factor and the Darcy-Weisbach friction factor. 

Equation: 


=->, Stanton — Pannell friction factor 


fe =?™, Fanning friction factor 
pw = =, Darcy — Weisbach friction factor (Moody) 


Um = (v), mean velocity 


The Reynolds number is expressed as a ratio of inertial stress and shear stress. 
Equation: 


pu, _— pumD — 2punR 


Nre = = 
HLUm/D id ld 


Both the friction factor and the Reynolds number have as a common factor, the kinetic energy 
per unit volume pu?, . This factor may be eliminated between the two equations to express 
the friction factor as a function of the Reynolds number. 

Equation: 


Tw 


_ 1 
fsp ~ Nrg bUm/D 


fr = 2 Tis 
Nre Um/D 
Tw 


—_ _8 
fow ~ Neg pUm/D 


Recall the expressions derived earlier for the wall shear stress and the average velocity for a 
Newtonian fluid and substitute into the above expressions. 


Equation: 
8 
f SP — Wr 
_ _16 
fr = "Nre 
64 
fow =~ Wre 


Correlation of friction factor versus Reynolds number appear in the literature with all three 
definitions of the friction factor and usually without a subscript to denote which definition is 
being used. 


Non-Newtonian fluids. The velocity profiles above were derived for a Newtonian fluid. A 
constitutive relation is necessary to determine the velocity profile and mean velocity for non- 
Newtonian fluids. We will consider the cases of a Bingham model fluid and a power-law or 
Ostwald-de Waele model fluid. The constitutive relations for these fluids are as follows. 
Equation: 


The Bingham Model 


dv 
Tee he a sete. ib Hag he 
dv, 
iy = 9 if eee | ee 
The Ostwald — de Waele (power — law) Model 
= dv, |” 1 dvuz 
Tap dy dy 


The power-law model is an empirical model that is often valid over an intermediate range of 
shear rates. At very low and very high shear rates limiting values of viscosity are approached. 


Assignment 8.1 


Flow in annular space between concentric cylinders as function of relative translation, 
rotation, potential gradient, flow or no-net flow. Assume incompressible, Newtonian fluid 
with small Reynolds number. The outer radius has zero velocity. Parameters: 


outer radius 

inner radius, maybe zero 

potential gradient, may be zero 

translation velocity of inner radius, may be zero 
rotational velocity of inner radius 


net flow rate, may be zero 


a. Express dimensionless velocity as a function of the dimensionless radius and 
dimensionless groups. Plot the following cases: 


Table of cases to plot 


Ri/Re OP/0dz Uz V6l1 


1 0 #0 #0 


2 0.5 #0 0 0 #0 
3 0.5 =0 #0 0 #0 
4 0.5 #0 #0 0 0 

5 0.5 #0 #0 #0 #0 


b. What is the net flow if the inner cylinder is translating and the pressure gradient is zero? 
c. What is the pressure gradient if the net flow is zero? Plot the velocity profile for this 
case. 


Assignment 8.2 
Capillary flow of power-law model fluid. Calculate the following for a power-law model fluid 
(see hint in BSL, 1960). 


a. Calculate and plot the velocity profile, normalized with respect to the mean velocity for 
n = 1,0.67, 0.5, and0.33. 

b. Derive an expression corresponding to Poiseuille law. 

c. Derive the same relation between friction factor and Reynolds number as for Newtonian 
flow by defining a modified Reynolds number for power-law fluids. 


Steady film flow down inclined plane 


Steady film flow down an inclined plane corresponds to case 4 (Churchill, 1988) or Section 
2.2; Flow of a Falling Film (BSL, 1960). These flows occur in chemical processing with 
falling film sulfonation reactors, evaporation and gas adsorption, and film-condensation heat 
transfer. It is assumed that the flow is steady and there is no dependence on distance in the 
plane of the surface due to entrance effects, side walls, or ripples. The Reynolds number must 
be small enough for ripples to be avoided. The configuration will be similar to that of BSL 
except x = 0 corresponds to the wall and the thickness is denoted by h rather than 6. 


It is assumed that the gas has negligible density compared to the liquid such that the pressure 
at the gas-liquid interface can be assumed to be constant. The potential gradient in the plane 
of the film is constant and can be expressed either in term of the angle from the vertical, 6 , or 
the angle from the horizontal, a. 

Equation: 


—-VP=pgcos B=pgsina, a=7/2-—8 


The equations of motion are as follows. 
Equation: 


_ @P 
C5 


dT zz 


0=pg cos B— = 


Uz 


0=pg cos B+ poe, Newtonian fluid 


Bo 


The boundary conditions are zero stress at the gas-liquid interface and no slip at the wall. 
Equation: 

v,=0, «r=0 

Fo =O;. 0h 


The shear stress profile can be determined by integration and application of the zero stress 
boundary condition. 
Equation: 


Tez = pgh cos B(#-1), 0<2<h 
Tw = —pgh cos B 


The velocity profile for a Newtonian fluid can be determined by a second integration and 
application of the no slip boundary condition. 
Equation: 


h? 2 2 
v, = babe | te (¢)’], 0<a<h 


_ pgh? cosB 
Uz, max — Qu 


The average velocity and volumetric flow rate can be determined by integration of the 
velocity profile over the film thickness. 


Equation: 
__ pgh? cosB 
(vz) ~~ 3u 
Q __ pgW h' cos 
7 age 


The film thickness, h, can be given in terms of the average velocity, the volume rate of flow, 
or the mass rate of flow per unit width of wall ("= ph (v-,)): 


Equation: 
p — ,| deed _ of _ 8H Q fw DP 
pg cos B pgW cos B p’g cos B 


Unsteady viscous flow 


Suddenly_accelerated plate. (BSL, 1960) A semi-infinite body of liquid with constant density 
and viscosity is bounded on one side by a flat surface ( the xz plane). Initially the fluid and 
solid surface is at rest; but at time ¢ = 0 the solid surface is set in motion in the positive x- 
direction with a velocity U. It is desired to know the velocity as a function of y and ¢t. The 
pressure is hydrostatic and the flow is assumed to be laminar. 


The only nonzero component of velocity is v; = vz (y, t). Thus the only non-zero equation 
of motion is as follows. 
Equation: 


The initial condition and boundary conditions are as follows. 
Equation: 


vz —0, t=0, y>O0 


v, =U, y=0,t>0 
vz=90, youw,t>d0 


If we normalize the velocity with respect to the boundary condition, we see that this is the 
same parabolic PDE and boundary condition as we solved with a similarity transformation. 
Thus the solution is 

Equation: 


V—e-= verte( 2) 
V4ut/p 


The presence of the ratio of viscosity and density, the kinematic viscosity, in the expression 
for the velocity implies that both viscous and inertial forces are operative. 


The velocity profiles for the wall at y = 0 suddenly set in motion is illustrated below. 


Wall set into motion; t= 0.01 0.03 0.1 0.3 1 


Developing Couette flow. The transient development to the steady-state Couette flow 
discussed earlier can now be easily derived. We will let the plane y = 0 be the surface with 
zero velocity and let the velocity be specified at y = L. The initial and boundary conditions 
are as follows. 

Equation: 


W=0, €=0, Oye DL 
i, =0, y=0, t20 
vz =U, vet, t>0 


or 
= 0, t=0, -L<y<L 
i. =—-U, =-[, t>0 


V,.= UC, y=L, t>0 


It should be apparent that the two formulations of the boundary conditions give the same 
solution. However, the latter gives a clue how one should obtain a solution. The solution is 
antisymmetric about y = 0 and the zero velocity condition is satisfied. A series of additional 
terms are needed to satisfy the boundary conditions at y= +. The solution is 

Equation: 


Uz = U Wig {erfe| Se =] — erfe| | \ 


Avt 


Developing Couette Flow; t= 0.01 0.03 0.1 0.3 1 


Asignment. 8.3 

Flow of a fluid with a suddenly applied constant wall stress. This problem is similar to that of 
flow (—oo < x < oo, y > 0, t > 0) near a wall (—co < x < co, y = 0) suddenly set into 
motion, except that the shear stress at the wall is constant rather than the velocity. Let the 
fluid be at rest before t = 0. At time t = 0 a constant force is applied to the fluid at the wall, 
so that the shear stress Ty, takes on a new constant value 7, at y = 0 fort > 0. 


a. Start with the continuity and Navier-Stokes equations and eliminate the terms that are 
identically zero. Differentiate the resulting equation with respect to distance from the 
wall and multiply by viscosity to derive an equation for the evolution of the shear stress. 
List all assumptions. 

b. Write the boundary and initial conditions for this equation. 

c. Solve for the time and distance dependence of the shear stress. Sketch the solution. 

d. Calculate the velocity profile from the solution of the shear stress. The following 
equation will be helpful. 

e. Suppose the fluid is not infinite but rather there is another wall at a distance 2h away 
from the original wall and it was also set into motion but in the opposite direction with 
the same wall shear stress. What are the stress and velocity profiles for the fluid between 
the two walls? Sketch and express as series solutions. 

f. What are the steady state stress and velocity profiles for the problem of part (e)? Sketch 
and express as analytical solutions. 


Multidimensional Laminar Flow 


Reading assignment 
Chapter 4 in BSL, Transport Phenomena 


Lubrication and Film Flow 


We already had two examples of flow in gaps that could be a thin film; 
Couette flow, and the steady, draining film. Here we will see that when the 
dimension of the gap or film thickness is small compared to other 
dimensions of the system, the Navier-Stokes equations simplify relatively 
and simple, classical solutions are possible. In Chapter 6, we saw that when 
the dimension of the gap or film in the x3 direction is small and the 
Reynolds number is small, the equations of motion reduce to the following. 
Equation: 


In the above equations, the subscript, 12, denote components in the plane of 
the gap or film. When the thickness is small enough, one wall of the gap or 
film can be treated as a plane even if it is curved with a radius of curvature 
that is large compared to the thickness. 


Lubrication flow with slider bearings. (Ockendon and Ockendon, 1995) 
Bearings function preventing contact between two moving surfaces by the 
flow of the lubrication fluid between the surfaces. The generic example of 
lubrication flow is illustrated with the slider bearing. 


y=H(x) 
6L 
0 x 
te ae 
U 1 
(a) Dimensional (b) Nondimensional 


A two-dimensional bearing is shown in which the plane of y = 0 moves 
with constant velocity U in the x-direction and the top of the bearing (the 
slider) is fixed. The variables are nondimensionalised with respect to U, the 
length L of the bearing, and a characteristic gap-width, hg, so that the 
position of the slider is given in the dimensionless variables. Again, 
referring back to Chapter 6, the dimensionless variables for this problem 
may be the following. 


Equation: 
*_ 2 *_ _Y _ H(z) 
c= 7) eae h(x) = ho 
kw *_ v L 
Ge PUSS ee 
p*- hip 


uLU 


Henceforth, the variables will be dimensionless with the * dropped. The 
boundary conditions are as follows. 


Equation: 
u=1, =0, y= 
Uu—vVU, =0, y= h(x) 
P= 0, o=0, 1 


The pressure can not suddenly equal the ambient pressure as assumed here 
because entrance and exit effects, but these will be neglected here. In 
reality, there may be a high pressure at the entrance of the bearing as a 
result of the liquid being scraped from the surface. The low pressure at the 


exit of the bearing may result in gas flowing in to equalize the pressure or 
Cavitation may occur. 


The dimensionless equations of motion and continuity equation are now as 
follows. 
Equation: 


(=25- =P =P) 
= dP eka 
O= ae + ay 


Ou Ov __ 
oT ap 


Integration of the equation of motion over the gap-thickness gives the 
velocity profile for a particular value of z. 
Equation: 


Notice that this profile is a combination of a profile due to forced flow 
(pressure gradient) and that due to induced flow (movement of wall). The 
velocity may pass through zero somewhere in the profile if the two 
contributions are in opposite directions. This is illustrated in the following 
figure. 


Fixed bcanm Fixed 
a 
_— | ~ 
=e SS 
ee = ———— S 
— —_— a 
— ————————— a 
———_ 
CE A nme vere ree oo DRAG sae 
! 
_—_—e ————> 


Velocity profile in slider bearing (Middleman, 
1998) 


Integration of the continuity equation over the thickness gives, 
Equation: 


= hf oa ) 
0 = f(t #) dy 
h h 
fo dy + v|o 


h 
a fo $e dy 


The latter integral can be expressed as follows. 


Equation: 
h h y=h 
Jo Se dy = A fp udy— ou 
d_ ph(z) 
mae ~ udy 
0 


Substituting the velocity profile into the above integral gives us the 
Reynolds equation for lubrication flow. 


Equation: 
d(hdP\) dh _y 
dx \ 6 dz dx 


Integration of this equation gives, 
Equation: 


h® dP 
—— =h+C 
6 dz . 


A second integration gives, 
Equation: 


where the constant of integration has to satisfy the boundary conditions, 
Equation: 


Io Ma)+C dn! = 0 


h?(z') 


Pressure profile is slider bearing (Middleman, 1998). 
Kj Ay. A> LAs. tan =(Ay— Ap) fi 


The role of the lubrication layer is to maintain a separation of the two 
surfaces in the presence of a load such that asperities (roughness) on the 
surfaces do not make contact. The load on the bearing is equal to the 
integral of the normal stress over the bearing surface. If the change in gap- 
thickness is small compared to the length, as assumed here, the normal 
stress is approximately equal to the pressure. If the gap-thickness is 


monotone decreasing, the pressure will be greater than ambient pressure 
inside the bearing. However, if the gap-thickness in monotone increasing, 
the pressure will be less than ambient in the bearing and it will have no load 
bearing capacity. If the gap-thickness is not monotone, then the pressure 
may be greater than ambient in some places and less than ambient in other 
places. If the bearing is designed to be load bearing, then a long section of 
decreasing gap-thickness and a short section of increasing thickness is 
desired. If the bearing is designed to be a scraper as piston rings then both 
sections of changing thickness will be short as to limit the amount of liquid 
passing through the gap. Gas entering the low-pressure region at the exit of 
the bearing surface prevents bearing surface contact from negative 
pressures. 


The analysis for the slider bearing can also be used to design an apparatus 
for depositing a uniform coating of a liquid on a substrate. 


Coating blade 


Coating flow (Middleman, 1998) 


Squeeze films. When two objects approach each other in a fluid their 
relative velocity is slowed as the resistance increases for the fluid leaving 
the gap. Here the relative velocity of the surface is in the normal direction 
rather than in the parallel direction as in the case of the slider bearing. This 
type of flow is important in the coalescence of emulsion droplets or foam 
bubbles. In the case of coalescence, hydrodynamics govern the dynamics of 


the approach of the surfaces to each other until the thinning is accelerated or 
retarded by surface forces (i.e., disjoining pressure). 


We will derive the classical Reynolds drainage of a liquid between two 
parallel disks of radius R approaching each other. The configuration of the 
system is shown below. 


U(t) 


Schematic of film drainage 
between parallel disks 


The system is symmetrical about its axis and the mid-plane. The thickness, 
h, is one-half of the distance between the disks and the velocity of each 
disk, U, is one-half of the approach velocity of the two disks. This 
nomenclature may be awkward but with this nomenclature, the solution also 
applies to the thinning of a liquid film between a solid surface and a gas 
bubble having zero shear stress at the interface. The velocity of approach of 
the disks may not be constant but rather the force pressing the disks 
together may be constant. Because of the symmetry, we will analyze the 
upper half-space with cylindrical polar coordinates. The system is 
axisymmetrical so the independent variables are r, z, t. The equations of 
motion and continuity equation in cylindrical polar coordinates are 
Equation: 


=-#+40(4) (Re) 
h 


0 “+0 
o= 2 + pF + 0(4)*+ O(Re) 
i 


The boundary conditions are 
Equation: 


The partial differential equations do not have an explicit dependence on 
time as time only enters thorough the boundary conditions. ‘Thus the 
variables will be made dimensionless with respect to the time dependent 
boundary conditions for the purpose of solving the PDE. 

Equation: 


kK oT. Kz 
es Sa ame 
* vr h *K__ Uz 
Up ps. U2 = Te 
His h3P 

pR2U 


The dimensionless equations and boundary conditions with the * dropped 
are now 
Equation: 


0=-2, =P=P(r) 
dP O70, 

0S Sar aa 

1 


or 

P=. pe] 
,=%=0, z=0 
= 0, Ve HL. Zed 


Integration of v, with respect to z in the equation of motion and applying 
the boundary conditions results in the velocity profile across the film 
thickness. 

Equation: 


Integration of the continuity equation over the film thickness gives, 
Equation: 


0= fo FEC) + Ss ) dz 


1 
— to 12 (rv,) dz+ Velo 


Tr 


The velocity profile across the thickness is substituted into the above 
equation and the integration preformed. 
Equation: 


Integration and application of the boundary conditions give 
Equation: 


P- 


(24 


| vo 


The pressure and radius can now be converted to dimensional variables so 
we can see the dependence of the parameters. 
Equation: 


rn Eh 


The pressure distribution has a maximum at the center of the disk and 
decreases to zero at the outer radius of the disk. The pressure integrated 
over the area of the disk gives the force required to bring the disks together, 
each disk with a velocity U, when each disk is a distance h from the 
midplane. 

Equation: 


3 ruR*U 
8 


R 
a (r) rdr 73 


This expression can be turned around to express the velocity of each disk 
approaching each other when a force F' is applied. 
Equation: 


This result is the classical Reynolds (1886) velocity for the thinning of two 
parallel disks. 


If the applied force is constant, the above equation can be integrated to 
determine the time it takes to thin down from some initial thickness, h;. 
Equation: 


1 1 #16 Ft 


he pA? 38 op Rt 


It the initial thickness is large but unknown, then it can be assumed to be 
infinity with only a small error in the time to thin down to a small thickness. 
An explicit expression for the time to thin from infinite thickness to a 
thickness h is 

Equation: 


_ 3 muR* 
16 Fh? 


From this expression, we see that it will take an infinite time to thin to zero 
thickness. In reality, as the film becomes very thin, surface forces 
(disjoining pressure) will become important in accelerating or retarding the 
rate of thinning. If the surfaces are solid surfaces, contact will be made at 
high points (roughness) and the contact stresses may limit the thinning. 


Transient Drainage of a Vertical Film Earlier we treated the steady film 
flow along an inclined plane. Here we will consider the transient drainage 
of a film that has zero flux at the upstream boundary. This corresponds to 
the transient behavior immediately after the flow if liquid is shut off in the 
problem of the steady flow along an incline plane. We will treat the wall as 
if it was vertical. If it is inclined from the vertical, the acceleration of 
gravity in the solution will need to be multiplied by the cosine of the angle 
from the vertical. It is assumed that the film is thin enough for the Reynolds 
number to be negligible and there are no ripples. Also, it is assumed that the 
thickness is large enough that surface forces (disjoining pressure) are 
negligible. The fluid in the film is assumed to be incompressible and 
Newtonian and the surrounding fluid is assumed to have zero density and 
viscosity. The surface tension and surface viscosity are neglected. The 
initial thickness of the film is assumed to be a constant value, h;. Let z be 


the coordinate in the direction of the film flow and z the direction 
perpendicular to the wall. The equations of motion and continuity equation 
for thin films, discussed in Chapter 6 have several terms that can be 
neglected. The resulting equations and boundary conditions follows. 
Equation: 
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The first equation states that there is zero potential gradient over the 
thickness of the film. Because the surrounding fluid has zero density and 
the surface tension is neglected, the pressure in the film is equal to that of 
the surrounding fluid. Thus, 

Equation: 


P=p— pgz=Po- pgz 
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The velocity profile across the thickness of the film can be determined by 
integrating the second equation. 


Equation: 
pa eght eo). 
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The flux or flow rate per unit width of the film can be determined by 
integrating the velocity profile over the thickness of the film. 
Equation: 
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h 
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The continuity equation can be integrated over the thickness of the film. 
Equation: 
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The derivative can be taken outside of the integral with the addition of 
another term that cancels the term in the previous equation. 
Equation: 
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Thus the differential equation for the film thickness is 
Equation: 


This is a first order, hyperbolic partial differential equation with constant 
initial and boundary conditions. Time and distance can be combined into a 
single similarity variable. The trajectories of constant values of the 
dependent variable can be calculated from the PDE. 


Equation: 
= _ Oh dh 
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Since the origin of all changes in thickness occur at the origin, the equation 
can be integrated as straight-line trajectories for each value of thickness 
between zero and the initial condition. 

Equation: 
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This is the classical solution for transient film drainage. Thick films initially 
drain very rapidly but the rate of drainage slows as the film thins. 


Notice that where the thickness has thinned below the initial thickness, the 
thickness is independent of the value of the initial thickness. Also, notice 
that the solution does not have a characteristic time, length, or thickness. 
This suggests that the thickness, time and distance are self-similar. In fact 
these variables can be combined into a single variable. 

Equation: 
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h < h; 


The thickness normalized with respect to the initial thickness can be 
expressed as a function of a single similarity variable or if a system length 
is specified, it can be expressed as a function of the dimensionless distance 
and time. 
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One may be interested in the volume of liquid that remains on a vertical 
wall of length L after the film is everywhere less than the initial thickness. 
This can be determined by integrating the film thickness profile over the 
length of the wall. 

Equation: 


fi ndz= 4/4, for t> Ht 


Ypgt ? pgh; 
1 4 L 
Jo hdz*= 4/5577, where 2*= 2/L 


This expression shows that the amount of liquid remaining on a vertical 
wall is inversely proportional to the square root of time. This solution is 
valid only after the film has everywhere thinned below the initial thickness. 


Assignment 9.1 
Transient drainage of vertical film. 


a. Combine the independent variables and parameters as a dimensionless 
similarity variable. Plot the normalized thickness as a function of the 
similarity variable. 

b. Suppose the length of the system is L. Plot the profiles of the 
normalized thickness as a function of dimensionless distance for 
different values of dimensionless time, (¢ = eps : 1: 20). 


Laminar Boundary Layer 


The flow of an ideal fluid (inviscid and incompressible) in two dimensions 
can be calculated for many configurations through the use of potential flow 
and complex variables. At a solid boundary, the ideal fluid has zero flux 
(the normal component of velocity is equal to that of the solid) but the 
tangential velocity may not be equal to that of the solid. Real fluids have a 
finite viscosity and (with only a few exceptions) the tangential component 
of velocity is equal to that of the solid, i.e., the boundary condition of no- 
slip. In many cases, the "external flow" sufficiently far from a solid body 
can be modeled as that of an ideal fluid and the effect of finite viscosity 
effects the flow only near the surface of the body (and downstream of the 
body). These situations can be treated by application of the boundary layer 
theory. 


The underlying assumption in the boundary layer theory is that there is a 
very thin layer near the body where the gradient of the tangential velocity is 
very large due to the action of viscosity and the no-slip boundary condition. 
Elsewhere the effect of viscosity is assumed to be unimportant and can be 
modeled as inviscid or potential flow. 


The continuity equation and equations of motion were specialized in 
Chapter 6 with the assumption that the boundary layer thickness,6 , is small 
compared to the characteristic dimension of the body, J, i.e., 6/I << 1. 
The equations (known as Prandtl's boundary layer equations) and boundary 
conditions are recalled here. 


Boundary layer flow 


along a wall 
(Schlichting, 1960) 


Equation: 
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The remaining terms in the equation of motion and continuity equation are 
of similar magnitude if 


Vy = OU 6/L), — p® = pus, 
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The assumption that the boundary layer thickness is small compared to the 
characteristic length of the body requires that the Reynolds number be large 
compared to unity if the terms in the equation of motion are to be of similar 
magnitude. If L is to represent the distance, x, from the leading edge of the 
boundary layer, the above relation describes the thickness of the boundary 
layer as a function of distance from the leading edge. 

Equation: 


Flow past flat plate at zero 
incidence (Schlichting, 1960) 


Another quantity that is of interest in boundary layer theory is the local drag 
coefficient due to the wall shear stress (some definitions differ by a factor 
of 2). The mean drag coefficient is the average value of this quantity over 
the surface of the body. 

Equation: 
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Laminar flow along flat plate 


The classical system for the study of laminar boundary layers is the flow of 
a fluid in uniform translation past a flat plate. The free stream velocity is 
constant and the pressure gradient is zero. The classical solution to this 
problem is the doctoral thesis of H. Blasius (1908). The equation of 
continuity is satisfied exactly by expressing the velocity as the curl of the 


stream function. Since the system does not have a characteristic length, a 
similarity transformation makes it possible to combine the two independent 


variables (x, y) into a single independent variable, 7 = vl ee The 


equations reduce to a quasilinear third order ordinary differential equation 
for the dimensionless stream function. The solution is given as a series 
solution. Its derivation is tedious and will not be discussed here. The reader 
is referred to Schlichting (1960) for details. 


An alternative approach is to keep the velocity components as the 
dependent variables and approximately satisfy the continuity equation by 
expressing the transverse velocity component in the equation of motion as 
follows. 

Equation: 
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dy 
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This is the approach taken in BSL. They express the solution as a cubic 
polynomial in 77. 
Equation: 


n= y/6, O(a) = 4. 64, / 42 
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Analogy _with wall suddenly set in motion 


The previously mentioned solutions may be accurate solutions to the 
boundary layer equations but do not offer much physical insight. Here we 
will derive the solution to the boundary layer flow by using the flow due to 
a wall suddenly set in motion discussed in Chapter 8. Since the wall extends 
to infinity, there is no dependence on x and the equations of motion, initial 
condition, and boundary condition are as follows. 

Equation: 


Equation: 
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The solution derived by a similarity transform in Chapter 7 is, 


Equation: 
= verte( —¥—) 
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The coordinates can be transformed such that the plate is stationary and the 
fluid is initially in uniform translation past the plate. The initial condition, 
boundary conditions, and solution are then as follows. 

Equation: 


Uz = Us, t=0, y>0 
Vz = 0, y=0, t>0 
vz =U, yorw,t>d 


Vz = Ux. erf (4) 


This exact solution describes the diffusion of momentum from the 

uniformly translating fluid to the stationary plate. It describes growth of the 

boundary layer as a function of the similarity variable, 7 = Ta . There is 
V 

no convection of momentum because there is no dependence on x and the 

transverse component of velocity is zero. 


Suppose now that the plate is not doubly infinite but only exists along the 
positive x axis and the flow is in the positive x direction.. Since there is 
now dependence on the z coordinate, boundary conditions depend on x and 
the convective terms in the equations of motion no longer vanish. Now 
consider the steady state flow past this semi-infinite plate. Assume that the 
flow is undisturbed until xz = 0. The differential equations are the boundary 
layer equations with zero velocity gradients. The equations and boundary 
conditions are as follows. 

Equation: 
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Recall that the convection terms are the convective derivative of the x 
component of momentum. Thus the equation of motion can be expressed as 
follows. 


Equation: 
Ove a Duy 
Du, __ - a 
Dt ~~ Ody? 


This equation can be used to describe the diffusion of momentum along a 
streamline which originated at z = 0 at t = 0. If the transverse velocity is 
zero, each streamline would be at a constant value of y. With steady flow, 
there is a one-to-one mapping between z and ¢ along each streamline. 
However, this mapping is not the same for all streamlines because different 
streamlines have different velocities. Beyond the boundary layer, the 
velocity is the free stream velocity and at the surface of the plate, y = 0, the 
velocity is zero. We will make the assumption that the mapping for the 


entire boundary layer can be approximated by using the average of the free 
stream velocity and the velocity at the plate. 
Equation: 


t(x) = CINE 
Qn 
Uo0 


Also, we assume for the first approximation that the transverse velocity is 
zero such that each streamline is at a constant value of y. The diffusion 
equation now transforms into a parabolic PDE in x and y. 

Equation: 
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This mapping of distance along the plate and time is substituted into the 
expression for the diffusion of momentum to a stationary plate that was 
introduced into a uniformly translating fluid at ¢ = 0. 

Equation: 
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vy) = U6. erf a 
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This first approximation neglects the transverse velocity and assumes that 
the time since passage of the front of the plate corresponds to the average 
velocity of the fluid in the boundary layer. Although this solution is quite 
close to the exact, Blasius solution, it does not satisfy the continuity 
equation. We now derive the second approximation by application of the 
continuity equation. 
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The limiting value of the transverse velocity with this approximation is, 


lim te . Ss 
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This limiting value of the transverse velocity differs from the Blasius 
solution only by a coefficient of 0.865. 


The transverse velocity results in convection of momentum away from the 
wall. If the convection velocity was constant, then its effect can easily be 
taken into account with the solution to the convection-diffusion equation. 
However, the transverse convection increases from zero at the wall to the 
limiting value in the free stream. Thus it makes more sense to use the 
average transverse velocity between the wall and at a point in the boundary 
layer for substitution into the solution of the convective-diffusion equation. 
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The second approximation will include the transverse convection by 
substituting the average transverse velocity into the convective-diffusion 
solution. 
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The first and second approximations are compared with the Blasius exact 
solution and quadratic and cubic approximation in the following figure. 
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The inclusion of the transverse convection made an insignificant 
improvement in the velocity profile. Thus it will be neglected in the 


following. The drag coefficient is calculated from the first approximation. 
Equation: 
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This drag coefficient differs from the Blasius solution only by the 
coefficient of 0.332 in the exact solution. 


The evolution of the boundary layer velocity profile with equal increments 
of distance from the leading edge of the plate is illustrated in the following 
figure. It is suggested that the student execute the plate.m file in the 
boundary directory to view the movie of the evolution of the velocity 
profile and to examine the equations used in the calculations. 
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Blasius solution for boundary layer flow past a flat plate 


The previous solutions were instructive in that they illustrated the 
correspondence with the diffusion of motion from a plate to the bulk fluid. 
However, the approximate solutions did not exactly satisfy either the 
continuity equation or the equations of motion. Blasius used the stream 
function to exactly satisfy continuity equation. The equations of motions 
were simplified by the boundary layer assumption that the thickness of the 
boundary layer is small compared to the distance from the leading edge of 
the plate. Also, the pressure gradient is zero for this case. 

Equation: 


Assume that the dimensionless velocity profile can be expressed as a 
function of a similarity variable. 
Equation: 


The approximate solution had the following form: 


Equation: 
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This suggests a similarity variable: 
Equation: 


The equation of motion is expressed in terms of the stream function. 
Equation: 
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Assume that the stream function is a function only of the similarity variable. 
Equation: 
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The dimensionless stream function is expressed as a function of only the 
similarity variable. 
Equation: 
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Substituting the above equations into the equation of motion and 
cancellation of two terms results in the following ordinary differential 
equation. 

Equation: 
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This is a third order ODE with two conditions at 7 —> oo and one condition 
at . It is convenient to solve it as a set of first order ODEs with initial 
conditions, two of which are specified and the third adjusted such as to 
satisfy the condition at infinity. 
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This set of ODEs can be solved numerically by one of the ODE solvers and 
the initial value of f” iterated until the boundary condition at infinity is 
matched. The code of this calculation is in the boundary directory as files, 
blasius.m, blasiusf.m, and balsiusd.dat. 


Assignment 9.2: Boundary layer flow past a wedge 


1. Derive the equations for boundary layer flow past a wedge. Use a 
factor of in the denominator of the similarity variable to be in keeping 
with contemporary textbooks. 

2. Use the code in the boundary directory of the CENG 501 website to 
solve the Flakner-Skan equation. 

3. Plot the velocity profiles as a function of the similarity variable for 
different angles of the wedge relative to the approaching free-stream 


velocity. Replace the parameter beta with the angle in degrees. 
A. Illustrate the boundary layer thickness by plotting contour lines of 10 
5. Plot the equally spaced streamlines for the same cases. 


Assignment 9.3: Flow in a wedge with zero shear at 0 = 0. 

Start from the continuity and Navier-Stokes equation and derive the 
equations for the flow field near a corner for flow in a wedge of fluid with 
no slip on one side and zero shear stress along 6 = 0. List all your 
assumptions. 


. Derive expressions for the stream function, velocity, and pressure. 

. For what distance from the corner is the solution valid? 

. What normal stress is required to keep the = 0 surface flat? 

. If the surface of zero shear stress can sustain only finite normal stress, 
in which way will the surface deform? Recognize that the no-slip 
surface can travel in either direction. 

5. After deriving the equations, view and plot the flow field for various 

angles using wedge.m file in the creeping directory of the CENG 501 

website. 
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Assignment 9.4: Rise of a spherical, inviscid bubble in a liquid. 

Start from the continuity and Navier-Stokes equation and derive the 
equations for the flow field of a spherical, inviscid bubble rising in a liquid 
by buoyancy. List all assumptions. 


1. Derive expressions for the stream function, velocity, and pressure. 
2. Derive the expression for the terminal rise velocity. How does it differ 
from the case with no slip? 


Flow with Free Surface 


We have already encountered free surfaces in systems such as the drainage of a 
liquid along a wall. In this case the free surface was a material surface and the 
boundary condition was that of continuity of pressure and shear stress. The same 
boundary conditions would be used for wind-driven waves on water and the shape 
of the vortex formed when water drains from a bathtub. The dimensionless 
numbers of importance are the Reynolds number Nr. = pU L/, Froude number 


2 
Np, = U?/gL, and gravity number Ng = 7 4 = — . These mentioned 


systems are of a macroscopic scale compared to surface forces and rheology and 
thus surface tension, surface elasticity, and surface viscosity were not significant. 
However, when the system dimensions become about 1 cm or less surface forces 
are no longer negligible and play an important role in the shape of the interface and 
in transport processes. The capillary number Nog = U/o and Bond number 
NBo = pg L?/c introduced in Chapter 6 become important dimensionless groups 
that quantify the ratio of viscous/capillary and gravity/capillary forces. As the 
dimensions decrease to about 1 mm we are in the range of capillary phenomena 
where surface tension and contact angles become important (e.g., the rise of a 
wetting liquid in a small capillary). As the dimensions decrease to 1 wm we are in 
the colloidal regime and not only is capillarity a dominant effect but also particles 
have spontaneous motion due to Brownian motion and thin films display optical 
interference as in the color of soap films. When the dimensions decrease to the 
range of 1 nm, it is necessary to include surface forces due to electrostatic, van der 
Waals, steric and hydrogen bonding effects to describe the thermodynamics and 
hydrodynamics of the fluid interfaces. At this scale the phases can no longer be 
assumed to be homogenous right up to the interface. The overlap of the 
inhomogeneous regions next to the interfaces results in forces that either attract or 
repel the interfaces. 


Boundary Conditions at a Fluid-Fluid Interface 


Analysis of macroscopic systems usually assume the fluid-fluid interface to simply 
be a surface of discontinuity in the density and viscosity of the bulk phases with no 
discontinuity in stress, (i.e., continuous pressure and shear stress across the 
surface). If there is no significant mass transfer, the surface is also a material 
surface and thus follows the motion of the fluid particles at the surface. 


Systems with a length scale about a centimeter or less and having fluid-fluid 
interfaces can no longer neglect the discontinuity in stress across the interface. A 
momentum balance across the interface is needed to describe the stress at the 


boundary. Also, if the system has surface-active components that affect the surface 
tension and/or surface viscosity, then a material balance is also needed to 
determine the composition of the interfacial region. A general treatment of the 
momentum balance at fluid-fluid interfaces is given in Chapter 10 of Aris, Vectors, 
Tensors and the Basic Equations of Fluid Mechanics. We summarize the results 
here assuming no slip at the surface and a Newtonian surface constitutive relation. 
The terms in the momentum balance are given on the left side and its description is 
given on the right side. 
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The surface constitutive equation for a Newtonian interface is (Slattery, 1990) 
TS = |o+ (kK +)V, e ve|IS + 2ee, 


The surface properties are a function of the composition of the interface. The 
species mass balance at the interface is given as (Edwards, et al., 1991) 
Equation: 


Or, 
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The parallel between the mass and momentum balances and constitutive equations 
for interfaces and bulk fluids should be noted (Gurmeet Singh, 1996). This analogy 
with bulk fluids is more complete if a surface pressure due to the reduction of the 
surface or interfacial tension due to adsorption is defined. 

Equation: 


The surface properties are as follows. 


Symbol Property 

Y mass per unit area 

Ln surface excess concentration of species n 
or surface or interfacial tension 
1 surface pressure 

K surface dilatational viscosity 
E surface shear viscosity 

a bag t surface tensors 

eG surface permutation symbol 
i. surface velocity 

J, surface diffusional flux 

A mean curvature 


K Gaussian or total curvature 


The normal component of stress. The discontinuity in the normal component of the 
total stress tensor for a hydrostatic system is as follows. 
Equation: 


[T], = pjn=2Hon 


where # is the mean curvature of the surface and a is the surface or interfacial 
tension. The mean curvature can be expressed as a function greatest and least 
curvatures of curves in the surface k, and kg or the principal curvatures in the 
directions of the principal curvature. 

Equation: 


2H = «1+ Ko 
The curvature in one of the principal directions can be expressed in terms of the 


equation for the arc of the curve. Suppose one principal direction is in the x-y 
plane. 
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Thus the curvature can be determined from the coordinates of the surface. We see 
that when the slope is small, the curvature can be approximated by the second 
derivative. A surface that is translational symmetric has zero curvature in one 
direction (e.g., surface of a cylinder). The two principal curvatures are equal on the 
surface of a sphere. A saddle shaped surface can have zero curvature because the 
principal curvatures have opposite signs. 


The difference in pressure across a curved interface is called the Laplace pressure 
after the Laplace-Young equation. This is a classical equation used to determine 


the shape of a static meniscus or to determine the surface tension from the shape of 
a Static meniscus such as the pendant drop shown here. This is a drop of water 
suspended from the tip of a capillary tube surrounded by air. Water surrounded by 
oil would be similar. Suppose we let the origin our coordinate system be the 
bottom of the drop. The system is an axisymmetric and has two radaii of curvature. 


Pendant water drop in air 
(Adamson, 1990) 


Equation: 
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The system is hydrostatic so the pressure profile can be expressed as the pressure 
jump at the origin and the difference in hydrostatic pressure along the profile. 
Equation: 


[p] = Ap. — Apgz 
Ap, = 20/r> 
lp] = 20/r. — Apgz 


where Ap, Laplace pressure at the apex of the drop and 1, is the radius of 
curvature of the apex of the drop. 


The surface or interfacial tension is found using the pendant drop analysis by 
estimating the value of the tension that best fits the calculated meniscus shape with 
the shape captured in a video image. 


When the pendent drop apparatus is designed so the meniscus is pinned to either 
the inner or outer edge of the tip of the needle, the contact angle does not influence 
the shape of the meniscus. If the drop rests on a flat surface, it is called a sessile 
drop and the elevation of the profile from the surface is a function of the contact 
angle that the meniscus makes with the substrate. 


A characteristic length scale can be determined for the hydrostatic meniscus. 
Suppose that the datum of elevation corresponds to the apex of the drop. The 
hydrostatic profile is described by the following dimensionless Young-Laplace 
equation. 

Equation: 
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If the dimensionless group is specified to equal unity, a characteristic length is 


defined. 
o 
L= 
V Apg 


Equation: 

This characteristic length, called the capillary constant, is sometimes defined with 
a factor of 2 multiplying the surface tension. It has a value of 2.7 mm for the water- 
air interface at ambient conditions. This length is representative of the meniscus 
height of water next to a vertical wall that is wetted by water. 


The dimensionless differential equations are now as follows. 


Equation: 
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The system of ODEs is computed from the apex of the drop where the 
dimensionless radius of curvature at the apex is a parameter, B. The value of B is 
adjusted until the best fit of the measured drop shape is obtained. 


Assignment 11.1 

Estimation of surface tension. Estimate the surface tension of a dilute surfactant 
solution, given the shape of a pendant drop shown here and data stored in the file, 
shape.dat in the Surface directory. You may use the code in pendant.m in the 
Surface directory to generate the dimensionless pendant drop profiles. Use units.m 
to plot the drop shape. Note: The program is not automated. Pendant.m computes 
the dimensionless shape and units.m converts it to cgs units to compare with the 
measured profile. 


Shape of Pendant Drop 


Surface tension gradients. If the system has only two components, i.e., the 
components comprising each phase then the surface or interfacial tension and 
contact angle is all that is required to describe the surface effects. However, if the 
system has another component that is surface active as to adsorb at the interface 
and reduce the surface or interfacial tension, then the interface must be treated as a 
two dimensional phase for which a mass and momentum balance is required. The 
mass per unit area of the surface-active component is the surface excess 
concentration or the amount adsorbed. If the system is not at equilibrium, i.e., not 
hydrostatic, then concentration gradients may exist in the interface that result in 
surface tension gradients in the interface. The difference between the clean 
interface tension and the local tension is called the surface pressure. The gradients 
in the surface pressure contribute to tangential stresses in the interface. It has the 
same role in the surface momentum balance as pressure gradient in the momentum 
balance for three-dimensional flow. For example, as liquid drains from a soap film, 
the drag of the liquid on the interface stretches the interface. The resulting 
expansion of the interface reduces the surface concentration of soap on the 
interfaces. This establishes a surface tension or surface pressure gradient between 
the interface in the film and the interface in the meniscus. This gradient tends to 
oppose the motion of the interface and thus tends to maintain the interface 
immobile as the liquid drains from between the two near-immobile interfaces. At 
the same time this surface tension gradient tends to pull the interface from the 
meniscus back into the film. This leads to a turbulent motion of the interface at the 
boundary between the meniscus and the film. This effect called "marginal 
regeneration" is a Marangoni effect caused by the surface tension gradient. 


Plate IV. An irregular 
mobile film in a vertical 
rectangle 


Plate I. Rigid film in a 
rectangular vertical 


Surface viscosity. Adsorption of a surface-active component at an interface not 
only changes the surface tension or surface pressure but can also affect the surface 
rheology. Material adsorbed at interfaces form two-dimensional surface phases that 
may be gasous, expanded liquid, condensed liquid, or solid. The surface viscosity 
can change by more than a order of magnitude at a transition from one surface 
phase to another. This is analogous to the change in viscosity of bulk fluids at 
phase transitions. The attached figure shows vertical soap film drainage of a 
system is similar to that of the mobile film except that dodecanol was added to the 
sodium dodecyl sulfate (SDS) solution. The dodecanol screens the electrostatic 


repulsion of the SDS at the interface and promote the formation of a condensed 
liquid monolayer. This monolayer is rigid in this system and the films drains much 
more slowly than in the case of the mobile film. The mechanism of this difference 
in the drainage of foam films has been explained in terms of the surface tension 
gradient driven instability and the stabilizing effect of a large surface viscosity 
(Joye, et al., 1994, 1996). 


Film Drainage and Deposition with Laplace Pressure 


In Chapter 8 we modeled the gravity drainage of a film along a wall neglecting the 
pressure between the liquid and gas because the mean curvature of the system was 
very small compared to the length of the film. Now suppose that we have a film 
that in connected to a curved meniscus. The meniscus may be moving along a 
substrate and depositing a film or gathering up a deposited film, e.g., a bubble ina 
capillary tube. Alternatively, the substrate may be stationary with respect to the 
meniscus and the film is draining into the mensicus, e.g. foam or emulsion film 
between two bubbles or drops. For simplicity, we will assume that we have a pure 
system so there are no surface tension gradients or effects of surface viscosity. We 
will assume that the system is translational invarient. A schematic of some possible 
system configurations are shown below. 


Section of a bubble in a horizontal tube. 


The transition region. 


Configration of a bubble in a tube (Breatherton, 1961) 


lubrication layer 


bao 


Film between bubbles and bubble against 
plate 


The continuity equation and equations of motion were specialized for lubrication 
and film flow in Chapter 6. The equations to O(h,/L) or O(h,/L)2 are as 
follows. 

Equation: 


O(hU12) dh _y 


0212 Ot 
2 
0= —ViaP + par, r3<h 
oP 
== x3 <h 


P=p-—Apgz 


The systems with a solid substrate will have the boundary condition of no-slip at 
the solid boundary and zero shear stress at the pure-fluid interface. In the case of 
two bubbles or drops coming into contact, the mid-plane is a plane of symmetry 
and has zero shear stress. It will be assumed the fluid interface is immobile in this 
latter case. The variable, h, is the half-film thickness in this case. Since this case 
has zero shear on one surface and no slip on the other surface, the solution will be 
derived as for the cases with the solid substrate. The boundary conditions are as 
follows. 

Equation: 


{ 0 at 2z3;=0, for stationary substrate or bubble 
V2 = 


+U at x3=0, for substrate moving relative to meniscus 


Ov» = 0 


Ors at 73 =h 


The pressure is uniform across the thickness of the film so the velocity profile can 
be determined by integrating the equation of motion over the film thickness and 
applying the boundary conditions. 

Equation: 


x 0 

Di (3 — has) ss fay 
0 

V12 = -£-VP+ a 


The average velocity is substituted into the equation of continuity. 
Equation: 


Ot - 3y 02x12 we U 


Oh 1: O(hPVi2P) | 0 
+ 


The pressure can be expressed in terms of the thickness by application of the 
Young-Laplace equation assuming that the system has no dependence on 3. 
Equation: 


o hit(x) 
sa + Segz 


a= / 
{1+ [r' (2)? 


Substituting the pressure into the previous equation results in a forth order 
quasilinear partial differential equation for the film thickness. These equations can 
be solved directly by numerical simulation. The equations will be specialized for 
special cases. 


Drainage of foam film. Deriving the equations in cylindrical coordinates and 
dropping the terms that originated from the substrate velocity and gravity can 
describe the drainage of the thin, horizontal, circular film between two bubbles. If 
the interface remains plane-parallel, the Reynolds film drainage model in Chapter 
9 can describe the drainage. The shape if the interface described by the above 
equations is "dimpled" and was investigated by Joye, et al., (1992). The pressure 
gradient in the film would be zero if the film were flat. The change in curvature 
where the film merges into the meniscus results in a large pressure gradient. This 
pressure gradient drains liquid from this region and this drainage causes a local 
thinning of the film. This leaves a thicker film or "dimple" in the center of the film. 
This axisymmetric drainage is unstable if the surface shear viscosity is small and 
the dimple will slip out to one side (Joye, et al., 1994, 1996). 
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1/2 thickness am 


0.00 0.05 0.10 0.15 


Film profile: R, = 1.8mm, 
Q = 107'°m?/8, y = 30dyn/cm, 
disjoining pressure not included. 


Bubble in a capillary. The motion of a bubble in a small capillary was investigated 
by Bretherton in 1961. With steady state, absence of gravity, and small slope in the 
film, the differential equations for the thickness of the film simplify to 

Equation: 


Bh  3uU h—b 


dx o h3 


where 0 is the asymptotic thickness of the film that is left on the wall far from the 
meniscus. Near the front of the film where the thickness is very large compared to 
the asymptotic thickness, the right-hand side of the above equation can be 
approximated by zero. The general solution in this region for some constants A, B, 
and C’ is as follows. 

Equation: 


Ph 0), h>>b 


dzx3 


he pa(mzy2 + B( mt)» +Cb 


This thick portion of the film merges with the spherical cap at the front of the 
bubble. Thus the asymptotic film thickness, 6, can be determined by requiring this 
portion of the film to have the same mean curvature as the front of the bubble. 
Equation: 


2/3 
9He thitea (22) 4 + +) near front of film 


dx? R a 
2H 4, in front of bubble 
thus 
2/3 
oe 3uU 
ke A(27) 


The constants of integration of the general solution of this film profile was 
determined matching with the shape of the meniscus at the front of the bubble with 
the solution by numerical integration. The asymptotic thickness of the film 
deposited by the advancing meniscus thus was found to be 

Equation: 


2/3 
DP: 0.043 ( 2) 
R oO 


where F is the mean radius of curvature of the spherical cap at the front of the 
bubble. 


In the thin film region, C'D, the film thickness is approximately equal to the 
asymptotic thickness and the differential equation can be linearized. 
Equation: 

Bh ~ 3nU h-b 


dx? ~— o b3 


with the general solution 


A =1+ae§+ Be cos (v3E/2) +e? sin (v3e/2) 
es 


The asymptotic film thickness is approached if the dimensionless length of the 
bubble is large compared to unity. It is assumed that this is the case. Different parts 
of the general solution apply to the front and back of the bubble. 

Equation: 


—l=ae, near front of bubble 


h 
b 
a —1= Be-? cos (v3¢/2) tye 7 sin (v3e/2), near rear of bubble 


The integration constant for the front part of the film, a , can be set to unity by a 
suitable choice of the origin of the axis. One integration constant for the rear of the 
bubble can be set to unity by suitable choice of the origin of the coordinate axis but 
the other must be specified to match the thickness of the film that is approaching 
the rear of the bubble. If this thickness is the asymptotic thickness, b, then there is 


a unique solution for the film at the back of the bubble. These profiles shows 
oscillations in thickness as the rear of the bubble is approached. If the approach 
thickness is b, then the minimum thickness is 0. 7166. The capillary suction at the 
back of the bubble results in a local thinning of the film similar to the local 
thinning that occurs in the foam film between two bubbles. 


EFFECTS IN SMOOTH CAPILLARY 
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Under static conditions, the Laplace pressures at the front and back of the bubble 
are both equal to 20/R with opposite sign, so the net pressure drop is zero. Under 
dynamic conditions the film profiles and thus the curvature at the front and back of 
the bubble are different. The dissimilar shapes of the front and back of the bubble 
results in a net dynamic pressure drop across the bubble. 


INDIVIDUAL LAMELLAE 
SMALL CAPILLARY, R<< f, 


Bee 


SMOOTH CAPILLARIES: CORRELATION OF MODEL 
WITH EXPERIMENTS 


Base (CALCULATED), cp 


Equation: 


Ap = 4.52( 2) (2) 


Oo 


Train of Bubbles or Foam in a Smooth Capillary. The resistance to flow of a single 
bubble in a tube is the starting point for the flow of a coarse foam through a tube. 
The contributing factors to the pressure drop are the slugs of liquid, interface 
deformation as described by Bretherton, and the additional resistance due to 
surface tension gradients. (Hirasaki and Lawson, 1985). Assume that the tube is 
small enough that the bubble train consists of bubbles separated by individual 
lamella. In this case the radius of curvature of the meniscus is smaller that the 
radius of the capillary. However, it can be determined from the capillary radius, 
bubble size, and foam quality or gas fraction. The comparison of the apparent 
viscosity for the flow of foam through the capillary tube is shown in the attached 
figure. These results show that interfacial deformation alone greatly underestimate 
the resistance to flow. Analysis shows that surface tension gradients have a 
significant effect in retarding the motion of the interface of bubbles when surface 
active material is present. The following figure shows the dimensionless interfacial 
velocity for different values of the dimensionless surface tension gradient group. 
The final figure shows that satisfactory agreement with theory can be obtained if 
surface tension gradient effects are taken into account. 
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Assignment 11.2 

Thickness of entrained film Compare the thickness of the film entrained by a plate 
pulled vertically out of a bath of liquid assuming one or the other of the two 
assumptions: (1) surface tension has no effect and (2) the capillary number is much 
less than unity. Make the thickness dimensionless with respect to a common 
reference length and compare the results as a function of the capillary number. You 
may accept the solution by Landau and Levich. 


Numerical Simulation 


The classical exact solutions of the Navier-Stokes equations are valuable for the insight they give and as an 
exact result that any approximate method should be able to duplicate to an acceptable precision. However, if 
we limit ourselves to exact solutions or even approximate perturbation or series solutions, a vast number of 
important engineering problems will be beyond reach. To meet this need, a number of computational fluid 
mechanics (CFM) codes have become commercially available. A student could use the CFM code as a virtual 
experiment to learn about fluid mechanics. While this may be satisfactory for someone who will become a 
design engineer, this is not sufficient for the researcher who may wish to solve a problem that has not been 
solved before. The user usually has to treat the CFM code as a "black box" and accept the result of the 
simulation as correct. The researcher must always be questioning if any computed result is valid and if not, 
what must be done to make it valid. Therefore I believe Ph.D. students should know how to develop their own 
numerical simulators rather than be just a user. 


We will start with a working FORTAN and MATLAB code for solving very simple, generic problems in 2-D. 
The student should be able to examine each part of the code and understand everything with the exception of 
the algorithm for solving the linear system of equations. Then the student will change boundary conditions, 
include transient and nonlinear capabilities, include curvilinear coordinates, and compute pressure and 
stresses. 


The Stream Function - Vorticity Method 


Two-dimensional flow of an incompressible, Newtonian fluid can be formulated with the stream function and 
vorticity as dependent variables. 


Equation: 
Ow Op _ a 
Ox? Oy? ay 
Dw _ dw Ow Ow __ O?w O?w 
Dt — a + Ye ae + Uy ay = (3 + ot | 
where 


V = (v2, ¥y 0) =VxA= (52, — 3, 0), A=(0,0,0) 


w= (0, 0, w) = (0, 0, 54 - $2) 


We will begin developing a generic code by assuming steady, creeping (zero Reynolds number) flow with 
specified velocity on the boundaries in a rectangular domain. The PDE are now, 
Equation: 


; O0<a<L,, 0<y<L,, Re=0 


This is a pair of linear, elliptic PDEs. The boundary conditions with specified velocity are 
Equation: 


Vnormal = specified 


Vitangential = Specified = O(U) 
vpco= f vends 


boundary 


WBC — Vx Missuaailg 


The boundary conditions at a plane of symmetry are 
Equation: 


nev=0 
neVv=0 
a = constant 
w=0 


The PDE and definitions are made dimensionless with respect to reference quantities such that the variables 
are of the order of unity. 


Equation: 
*__ _& *_ Ye Le 
eS ge! Ose ge! A Sh 
= _ »¥ = 
vey, ed, vt 2 
Yo_| O* Yo | Oye * 
Es Ox*? am En Oy*? a 
: Veo Cag 1 
Qay* Qay* 
oe ar a? ae =0 
p= wv | ay* v,* wo | av* 
x UL, | OyF? “¥ — UL, | Ox* 


We have four unspecified dimensionless groups and two unspecified reference quantities. We can specify two 
of the groups to equal unity and the other two are equal to the aspect ratio or its square. 


Equation: 
Do eae Yo — 
[as|=4 [ve] = 


=> Yo = L,U, Wo = 


oS 
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The PDE and definitions with the * dropped are now as follows. 


Equation: 
oO? o? 
= +a? =e 
; 0<a< 1, 0< y< 1 
(cha oO? 
Oa? + a? Oye =0 
Ow Ow 
Vz = ABy > Vy => — Bx 
__ Ovy Ov, 


WwW Be Ody 


Finite Difference Approximation 


The PDE and BC will be solved using a finite difference method. A grid point formulation rather than a grid 
block formulation will be used since the dependent variables are specified at the boundaries rather than their 
flux or normal derivative. The computation domain will be discretized such that the boundary conditions and 


dependent variables are evaluated at 21, %2,...,®7mi, CIMAX and Y1, Y2,---, YIM1, YIMAX- The finite 
difference grid appears as follows. 


Finite difference grid for discretizing PDE (grid 
point formulation) 


The unit square is now discretized into JMAX by JMAX grid points where the boundary conditions and 
dependent variables will be evaluated. The grid spacing is 6 = 1/(JMAX — 1). The first and last row and 
column are boundary values. 


The partial derivatives will be approximated by finite differences. For example, the second derivative of 
vorticity is discretized by a Taylor's series. 
Equation: 


_ (a & ( & 
Wi41 = w; + 6(S2), + $ (St) 
: = — §( ae & (aw) _ & ( Fw 
Wi-1 = Wi 6(#),+ 2 (27), 6 (27). 
WwW 


aw _ wi-1—2 Wi tw _ oo 
Ox? J, & 24 


= int Oi + O(8?) 


The finite difference approximation to the PDE at the interior points results in the following set of equations. 
Equation: 

bisig + Wi-rg + rdiju + ar piz-1 — 2 (1 + a?) wi + &wi; =0 

Wiig + wi-1g + 7? wij + OP wi j-1 —2(1+a7)w,; =0 

i,j =2, 3, ---JMAX-1 


The vorticity at the boundary is discretized and expressed in terms of the components of velocity at the 
boundary, the stream function values on the boundary and a stream function value in the interior grid. (A 
greater accuracy is possible by using two interior points.) The stream function at the first interior point 
(¢ = 2) from the x boundary is written with a Taylor's series as follows. 

Equation: 


ve = w£0(3b) +4 (Sr), +0 (8) 
2 
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The choice of sign depends on whether z is increasing or decreasing at the boundary. Similarly, on a y 
boundary, 
Equation: 


2a? (2 — ve?) ZaveO — AvBC 
BC __ 1 £ y 
<= 52 ca ag tO) 


The boundary condition on the stream function is specified by the normal component of velocity at the 
boundaries. Since we have assumed zero normal component of velocity, the stream function is a constant on 
the boundary, which we specify to be zero. 


The stream function at the boundary is calculated from the normal component of velocity by numerical 
integration using the trapezodial rule, e.g., 


bj = 05-1 + £[(v2);1 + (v2),] /0 


Solution of Linear Equations 


The finite difference equations for the PDE and the boundary conditions are a linear system of equations with 
two dependent variables. The dependent variables at a x;, y; grid point will be represented as a two 
component vector of dependent variables, 

Equation: 


The pair of equations for each grid point can be represented in the following form 
Equation: 
aij Wisi + bij Wi-1y + C1 Wig + dig Wig + 4, Wig = fi; 
i,j = 2, 3,---JMAX-1 


Each coefficient is a 2x2 matrix. For example, 
Equation: 


E4511 Pig + €43,1,2 Wij 


e; i Uj7 = 
ae ap €1,5,2,1 Wig + €i,5,2,2 Wi, 


The components of the 2x2 coefficient matrix are the coefficients from the difference equations. The first row 
is coefficients for the stream function equation and the second row is coefficients for the vorticity equation. 
The first column is coefficients for the stream function variable and the second column is the coefficients for 
the vorticity variable. For example, at interior points not affected by the boundary conditions, 

Equation: 


The coefficients for the interior grid points adjacent to a boundary are modified as a result of substitution the 
boundary value of stream function or the linear equation for the boundary vorticity into the difference 
equations. The stream function equation is coupled to the vorticity with the e;; coefficient and the vorticity 
equation is coupled to the stream function through the boundary conditions. For example, ata z = 0 
boundary, the coefficients will be modified as follows. 

Equation: 


¢= 2 
fG,3.1) = f@G4,1)=0G, L077 

BC 
f (i, 5,2) = F (45,2) — B (i, 4,2, 2)* (2*UBC/5 — 2*vBC (3) /5 — 0S ) 
e (i, j,2,1) =e (i, 9,2, 1) — b (i, 5, 2,2)*2/e 


If the boundary condition is that of zero vorticity as on a surface of symmetry, it is necessary to only set to 
zero the coefficient that multiplies the boundary value of vorticity, i.e., 


i=2 
b(i, j,2,2) =0, when w(1,7) =0 


The linear system of equations and unknowns can be written in matrix form by ordering the equations with 


the i index changing the fastest, i.e., 7 = 2,3,..., JMAX — 1. The linear system of equations can now be 
expressed in matrix notation as follows. 
Equation: 
Aex=F 
where 
W2,2 
W2,2 
u22 
3,2 
3,2 
x = = W3,2 
UujM1,JM1 
WIM1,JM1 
WJM1,JM1 


The non-zero coefficients of coefficient matrix A for JM AX = 5 is illustrated below. Since the first and last 
grid points of each row and column of grid points are boundary conditions, the number of grid with pairs of 
unknowns is only 3x3. Only one grid point is isolated from boundaries. The 7,7 location of the coefficients 
become clearer if a box is drawn to enclose each 222 coefficient. 


nz = 83 


Most of the elements of the coefficient matrix are zero. In fact, the coefficient matrix is a block pentadiagonal 
sparse matrix and only the nonzero coefficients need to be stored and processed for solving the linear system 
of equations. The non-zero coefficients of the coefficient matrix are stored with the row-indexed sparse 
storage mode and the linear system of equations is solved by preconditioned biconjugate gradient method 
(Numerical Recipes, 1992). 


The row-indexed sparse matrix mode requires storing the nonzero coefficients in the array, sa(k), and the 
column number of the coefficient in the coefficient matrix in the array, ija(k), where k = 1,2,.... The 
indices that are needed to identify the coefficients are the z, 7 grid location (2, 7 = 2,3,..., JM AX — 1), the 
m equation and dependent variable identifier (m = 1 for stream function; = 2 for vorticity), and the [A row 
index of the coefficient matrix (IA = 1,2,3,4,...,2(JMAX — 2)2). The coefficient matrix has two 


equations for each grid point, the first for the stream function and the second for vorticity. The total number of 
equations is NN = 2(JMAX — 2)2. 


The diagonal elements of the coefficient array are first stored in the sa(k) array. The first element of ija(k), 
ija(1) = NN + 2 and can be used to determine the size of the coefficient matrix. The algorithm then cycles 
over the pairs of rows in the coefficient matrix while keeping track of the 2,7 locations of the grid point and 
cycling over equation m = 1 and 2. The off-diagonal coefficients are stored in sa(k) and the column 
number in the coefficient matrix stored in ija(k). As each row is completed, the & index for the next row is 
stored in the first NN elements of 77a. 


A test problem for this code is a square box that has one side sliding as to impart a unit tangential component 
of velocity along this side. All other walls are stationary. The stream function, vorticity, and velocity contours 
for this problem with a 20x20 grid is illustrated below. 


Stream function contours for square 
with top side sliding 


Velocity Arrays 


Velocity vectors for square with top 
side sliding 


Assignment 12.1 

Use the code in the /numer/ directory of the CENG501 web site as a starting point. Add the capability to 
include arbitrarily specified normal component of velocity at the boundaries. The stream function at the 
boundary should be numerically calculated from the specified velocity at the boundary. Test the code with 
Couette and Plane-Poiseuille flows first in the z and then in the y directions. Display vector plot of velocity 
and contour plots of stream function and vorticity. Show your analysis and code. 


Assignment 12.2 
Form teams of two or three and do one of the following additions to the code. Find suitable problems to test 
code. 


1. Calculate pressure and stress. Test problems: Couette and Plane-Poiseuille flows. 

2. Calculate transient and finite Reynolds number flows. Test problems: Plate suddenly set in motion and 
flow in cavity. 

3. Add option for cylindrical polar coordinates. Use the coordinate transformation, z =In r . Test problems: 
line source and annular rotational Couette flow. 


Form new teams of two or three with one member that developed each part of the code. Include all features 
into the code. Test the combined features. Work together as a team to do the following simulation 
assignments. 


Assignment 12.3 
Boundary layer on flat plate. Compute boundary layer velocity profiles and drag coefficient and compare with 
the Blasius solution. What assumption is made about the value of the Reynolds in the Blasius solution? 


Assignment 12.4 
Flow around cylinder. Assume potential flow far from the cylinder. Calculate drag coefficient as a function of 
Reynolds number and compare with literature. 


Assignment 12.5 
Flow around cylinder. Find the Reynolds number at which boundary layer separation occurs. 


FORTRAN on Owlnet 


Modern FORTRAN compilers in the MS Windows are very user friendly. They have built in project 
workspace routines, documentation, and integrated editor and debugger. If you use the FORTRAN compiler 
on Owlnet, it is useful to know how to use make files to compile and link FORTRAN code. The following is 
an example make file. You may give it the name, makefile, with no extension. 


f77 -c streamfi.f; 
f77 -c be1.f; 

f77 -c coefi.f; 
f77 -c asolve.f; 
f77 -c atimes.f; 
f77 -c dsprsax.f; 
f77 -c dsprstx.f; 
f77 -c linbcg.f; 
f77 -c snrm.f; 

f77 - 

o exe streamf1.o bc1.0 coef1.o asolve.o atimes.o dsprsax.o dsprstx.o linb 
cg.oO snrm.o 


You have to give yourself permission to execute the make file with the one-time command, 
chmod +x makefile 


The line with the -c option compiles the source code and makes an object code file with the extension of .o. 
The line with the -o option links the object files into the executable file called exe. This executable file can be 
executed by typing its name, exe, on the command line. 


This example makefile is not very efficient because it will compile every source code listed. There are 
instructions to recompile only the source code that has been modified or "touched" since the last time the 
makefile was invoked. However, I do not recall the instructions. 


If you are not familiar with FORTRAN, you should get a paperback book on FORTRAN-77 or FORTRAN- 
90. An example is D. M. Etter, Structured FORTRAN 77 for Engineers and Scientists. The following is a 
tutorial from the CAAM 211 webpage. http://Awww.caam.rice.edu/ caam211/JZoss/project1. html 


Computer facilities that have compilers usually have an interactive debugger. However, I have not found how 
to access the debugger on owlnet. The debugger with the Microsoft FORTRAN Powerstation (now HP 
FORTRAN) is very easy to use. 


Transients and Finite Reynolds Number 


The development of numerical simulators should start with a problem for which exact solutions exist so that 
the numerical algorithm can be verified. Steady, creeping flow has many exact solutions and thus it was used 
as the starting point for developing a numerical simulator for the Navier-Stokes equations. The student should 
verify that their code is able to duplicate exact solutions to any desired degree of accuracy. 


The next stage in the development of the code is to add transients and finite Reynolds numbers. Finite 
Reynolds number results in nonlinear terms in the Navier-Stokes equation. Algorithms for solving linear 
systems of equations such as the conjugate gradient methods solves linear systems in one call to the routine. 
Nonlinear terms need to be iterated or lagged during the transient solution. Thus nonlinear steady state flow 
may be solved as the evolution of a transient solution from initial conditions to steady state. 


The transient and nonlinear terms now need to be included when making the vorticity equation dimensionless. 


Equation: 
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Two of the dimensionless groups can be set to unity by expressing the reference time in terms of the ratio of 
the characteristic velocity and characteristic length. The remaining dimensionless group is the Reynolds 
number. 

Equation: 


The dimensionless equation for vorticity with the * dropped is 
Equation: 


+ 


Ow Ow Ow » O?w 
oe = a 
Ox? Oy? 


The convective derivative can be written in conservative form by use of the continuity equation. 
Equation: 


(vj w;) ; = Vig Wj + Vi Wii 

but 

vi; =0, for incompressible flow 
thus 

05 W34 = (Vi W;) ; 
or 


vVw=Vevw 


We will first illustrate how to compute the transient, linear problem before tackling the nonlinear terms. 
Stability of the transient finite difference equations is greatly improved if the spatial differences are evaluated 
at the new, n + 1, time level while the time derivative is evaluated with a backward-difference method. 
Equation: 


writ? Awi; 
© a + 0 (At) = ae + 0 (At) 
tafe fo n+1 n 
Aw; Wij Wi 5 


Recall that the finite difference form of the steady state equations were expressed as a system of equations 
with coefficients, a, b, , ect. The vorticity equation for the 7,7 grid point will be rewritten but now with the 
transient term included. It is convenient to solve for Aw rather than w”*? at each time step. 


Equation: 


2 
aAw? + baw), +cAwtt! +dAuwti!, +e Aw! — SR Aw?tt 


+1) J tj+ tJ aj iJ 
sate fe n _ n = n _ n = n 
= awe — bw 5 — CWP — awe, — ew, 


We can see that this equation is of the same form as before with only the coefficients e and f of the vorticity 
equation modified. 
Equation: 


1. &Re 
e =e At 


— = n = nm — n Pe n a 1) 
fl=f aw; bw); C Wi j41 dw, 4 ew; ; 


Since the variable in the linear system of equations is now Aw rather than w, the equations for the stream 
function also need to be modified. 
Equation: 


ai Wisi + O11 Wi-1g + 11 Wager + dia Vijri + e11 Pig + e12 Awi; 
= fia —e12 wi; 
thus 


f' = fii —e12 we; 


These modifications to the coefficients should be made after the coefficients are updated for convection and 
boundary conditions. 


The numerical calculations will start from some initial value of vorticity and proceed to steady state. The size 
of the time step is important if accuracy of the transient solution is of interest. The truncation error of the time 
finite difference expression is approximately the product At — The time step size can be chosen as to keep 
this value approximately constant. Thus the time step size will be chosen as to limit the maximum change in 
the magnitude of w over a time step. Let Aw, be the maximum change in w over the previous time step 
and Aw spec be the specified value of the desired maximum change in w. The new time step can be estimated 
from the following expression. 


Equation: 


Aw spec 


Atnew = Abel Fion.. 


The initial time step size needs to be specified, and the new time step may be averaged with the old or 
constrained to increase by not greater than some factor. 


Now that we have a means of calculating the transients solution, this approach can be used to treat the 
nonlinear terms. Recall that when the equations were linear, the system of equations had constant coefficients 
and the conjugate-gradient linear solver solved the linear system of equations. 

Equation: 


Ax=F 


The convective derivatives have a product of velocity and vorticity, which can be expressed as a product of 
the derivative of stream function and vorticity. One approximation would be to use the value of the stream 


function from the old time step to calculate the coefficients for the new time step. An alternative is use a 
predictor-corrector approach, similar to the Runge-Kutta solution for quasilinear ordinary differential 
equations. The "predictor" step estimates the solution at the new time step using the coefficients calculated 
from the stream function of the previous time step. This gives an estimate of the stream function at the new 
time step. This estimated stream function is then used to evaluate the coefficients for the "corrector" step. The 
steps in the calculations are as follows. 

Equation: 


A (x”) xrtl* _ Fp (x”) 


A be) xvtl_F- Gere) 


The choice of the finite difference expression for the convective derivative is very important for the stability 

of the solution. We will illustrate here treatment of only the x derivative. Suppose 7 — 1,7, andz + 1 are the 

grid points where the stream function and vorticity are evaluated and i — 1/2 andi + 1/2 are the mid-points 
where the velocity-vorticity products are evaluated. 

Equation: 


Ouw (uw)isi/2 _ (uw);1/2 
dz Az 


The dependent variables are known only at the grid points and interpolation will be needed to evaluate in 
between location. If the product is evaluated at i — 1/2 andi + 1/2 by using the arithmetic average of the 
values at the grid points on either side, the numerical solution will tend to oscillate if the convective transport 
dominates the diffusive (viscous) transport. It is preferable to determine the upstream direction from the sign 
of Uj 41/2 OF U4_-1/2 and use the vorticity from the upstream grid point. This is illustrated below. 


ui-1/2 = 0.5*(u(i,j)+u(i-1,j)) 
if ui-1/2 >0 then 
(uw)i-1/2 = ui-1/2 wi-1 ; 
b'=b + $\delta$ Re ui-1/2 
else 
(uw)i-1/2 = ui-1/2 wi ; e'=e + $\delta$Re ui-1/2 
endif 


uiti/2 = 0.5*(u(i, j)+u(it1, j)) 
if uI+1/2 >0 then 
(uw)it+1/2 = uit1/2 wi ; 
e'=e - $\delta$ Re uit+1/2 
else 
(uw)it+1/2 = uit1/2 witi ; 
a'=a - $\delta$ Re uit+1/2 
endif 


The y-direction will be similar except the aspect ratio, , must be included. These modifications to the 
coefficients must be made before the coefficients are updated for the boundary conditions. 


Calculation of Pressure 


The vorticity-stream function method does not require calculation of pressure to determine the flow field. 
However, if the force or drag on a body or a conduit is of interest, the pressure must be computed to 
determine the stress field. Here we will derive the Poisson equation for pressure and determine the boundary 


conditions using the equations of motion. If the pressure is desired only at the boundary, it may be possible to 
integrate the pressure gradient determined from the equations of motion. The dimensionless equations of 
motion for incompressible flow of a Newtonian fluid are as follows. 


Equation: 
on ok * 
ae = oe + vteVt¥vt= oe + V*e(v*v*) = 
—V*Pt+ 3 Vv 
where 
ve= G, =F, te ft, Pach, V=LV 


P = p-—pgz, Re= 


Here all coordinates were scaled with respect to the same length scale. The aspect ratio must be included in 
the final equation if the coordinates are scaled with respect to different length scales. We will drop the * and 
use x and y for the coordinates and u and v as the components of velocity. 


The following derivation follows that of Hoffmann and Chiang (1993). The equations of motion in 2-D in 
conservative form is as follows. 


Equation: 
du , Aw) | dur) __ _ aP | 1 (du, du 
ot a Ox au Oy Ox + Re \ dx? F Oy? 
ov Our) a(v?) _ _ OP 1 haa 0?u 
ot ca Ox + oy Oy a Re \ dx? ce Oy? 


The Laplacian of pressure is determined by taking the divergence of the equations of motion. We will carry 
out the derivation step wise by first taking the x-derivative of the x-component of the equations of motion. 


Equation: 


6) (=) Ou Ou 07u —s- Ou Ov 02v Ov Ou 07u 0?2P 1 0 


i 2 = Vv? 
Ot \ Ox Ox Ox 7 nee 7 Ox iy Oxbu = Ox Oy 1 Baer Ox? TiRe aa | 4) 


Two pair of terms cancels because of the equation of continuity for incompressible flow. 
Equation: 
Ou ( Ou dv \) _ 
De (# + gu) = 0 
O?y Ou _,, af du dv \ _ 
U Brdy + Wo = US (2 ae gr) 0 
Thus 


a (du du \2 au , dv du au _ ss P 1 9/2 
a (oe) (oe) TU ger oh og ay! oxay = Ox? Re da (V"u) 


Similarly, the y-component of the equations of motion become 
Equation: 


0 (av dv)? Ou | Ou dv | Pu _ ae, 
Ot \ dy Oy Oy? = Oy Ox dxdy —S Oy? Re Oy 


The x and y-components of the above equations are now added together and several pairs of terms cancels 
pair-wise. 
Equation: 


a) 
Of Bu O*u oO f[ dv a’u\ _ 0? f du Ov a? [ du Ov \ _ 
Ox (33 cs ot | a & ($3 iF oH) ~~ Ox? (# a | + Oy? (% + | =0 
Therefore the equations reduce to 
Equation: 
Ou ae Ov "4 ou Ov _ O°P , oP 
Ox Oy Oy Ox Ox? Oy? 


The left-hand side can be further reduced by consideration of the continuity equation as follows. 
Equation: 


2 2 
Ou Ov du \2 Ov du Ov 
(3 + 3) = (a) + (#) + 235 ay = 9 


from which 


2 
du \2 Ov — Ou Ov 
(sr) a (#) _ 255 Oy 


Upon substitution, we have 
Equation: 


(SF +) -o(F Ov Ou =) 
Ox? dy?) ~\ dy Ox Ox Oy 


This equation can be written in terms of the stream function. 
Equation: 


OP oP _ [ab ay (o> \’ 
Ox? dy? | Ox? Oy? OxOy 


This is the equation when the coordinates have the same reference length. If the coordinates are normalized 
with respect to different lengths, then the equation is as follows. 
Equation: 


OP 10°P _,2| dy a (= ) 
Ox? Oy? Ox? Oy? OxOy 


The Poisson equation for pressure needs to have boundary conditions for solution. The equations of motion 
give an expression for the pressure gradient. The normal derivative of pressure at the boundary can be 
determined for the Neumann-type boundary condition. The equations of motion usually simplify at 
boundaries if the boundary has no slip BC or if it has parallel flow in the direction either parallel or normal to 
the boundary. For example, for no-slip 
Equation: 
2 2. 

oF =e ane at y=c, u=v=0 

oP 1 d?v 


Oy = Re Ba?” at r=c, u=v=0 


If the flow is parallel and in the direction normal to the boundary 
Equation: 


OP _ ao & 
Ge = Re Op? y Cc, U 0, u u(y) 


oP _ 1 
5, = Re Ba?? g=60=—0.0=0 (2) 


Cylindrical-Polar Coordinates 


The code for the numerical solution of the Navier-Stokes equations in Cartesian coordinates can be easily 
modified to cylindrical-polar coordinates. The coordinate transformation is first illustrated for the Laplacian 
operator. 

Equation: 


1 1 0? 
vba io (Se) + eee m1 <7Tr<ra, 0<0<4@ 


The independent variables are made dimensionless with respect to the boundary parameters. 
Equation: 


Ps). oS Ley <p ee et 
r 


1 


The Lapacian operator with the dimensionless coordinates after dropping the * is now, 
Equation: 


vp= aloo (re) CS 1l<r<f8, 0<0<1 


Lr or\" Or) 62 1? 00 


The radial coordinate is transformed to the logarithm of the radius. 


Equation: 
2= 25 =ylnr, 0<z< 1, 1 = inn) 
r =exp (z/7) 


~~ 2 @z2 


2 92 
LB (rg) = 4 ary 


The Laplacian operator is now as follows, 
Equation: 
1 [y? 07 a? O*y 1 


, O<z<1, 0<6<1, a=— 
ro lr? O22 | 2 0 < ane 


Vp = 


The finite difference expression for the Laplacian will be same as that for Cartesian coordinates with (z, @) 
substituted for (x, y) except for y? and r? factors. 
Equation: 


The curl operator is modified from that in Cartesian coordinates. 


Equation: 
—_ a wy 
Ur = “7 80 
= Oy _ oy OW 
V6 = or OF Oe 
1 O(r v9) a Ov, 
WS or r 06 
_ x Ar v9) a Ov, 
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The vorticity boundary conditions with the transformed coordinates is for the z boundary, 
Equation: 


2 BC 
Bo_ _ 7 2 __ ,1,BC y 2 BC a Ov; 
wp = “a gr (Y2 0 te a 
and for the 6 boundary, 
Equation: 
BC 
we = Sal — pBC) 4 a2 ec 7 O(rvBC) 
1 2 62 1 os r 5 r po az 


The stream function at the boundaries are expressed different from that in Cartesian coordinates. 
Equation: 

dy =(rv,/a)d6, at z boundary 

dy = —(r vg/y) dz, at 0 boundary 


The convective terms are expressed different from that in Cartesian coordinates. 
Equation: 


O(uw) 1 Aruw) 
Oz = r Or 
y Aruw) 
~~ 72 Oz 
ey x 
PF ea, (rww) 4/2 — (ruw);_1| 
O(vw) a vw) 
a Oy > r 00 
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The code should be written so one can have either Cartesian coordinates or transformed cylindrical-polar 
coordinates. A parameter will be needed to identify the choice of coordinates, e.g. icase = 1 for Cartesian 
coordinates and tcase = 2 for transformed cylindrical-polar coordinates. Also, another parameter should be 
specified to identify the choice of boundary conditions, e.g., 2bc = 1 for radial flow, 7bc = 2 for Couette flow, 
and ibe = 3 for flow around a cylinder. Test cases with known solutions should be used to verify the code. 
The first case is radial, potential flow from a line source and the second is Couette flow in the annular region 
between two cylindrical surfaces. 


Equation: 
Role. , 
radial flow 
ve = 0 
v, =0 


ve = (r —1/r)/(B- rh Couette flow 


Flow around a cylinder needs a boundary condition for the flow far away from the cylinder. The flow very far 
away may be uniform translation. However, this condition may be so far away that it may result in loss of 
resolution near the cylinder. Another boundary condition that may be specified beyond the region of influence 
of the boundary layer is to use the potential flow past a cylinder. This boundary condition will not be correct 
in the wake of cylinder where the flow is disturbed by the convected boundary layer buts its influence may be 
minimized if the outer boundary is far enough. 

Equation: 


v, = (1—1/r?) cos 6 
; Potential flow past cylinder 
vg =—(1+1/r?) sind 


Potential flow will not be a valid approximation along the @ boundaries close to the cylinder. Here, one may 
assume a surface of symmetry, at least for the upstream side. 


